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ABSTRACT 


Signal  processing  methods  for  signals  sampled  at  different  rates  are  inves¬ 
tigated  and  applied  to  the  problem  of  signal  and  image  reconstruction  or  super- 
resolution  reconstruction.  The  problem  is  approached  from  the  viewpoint  of  linear 
mean-square  estimation  theory  and  multirate  signal  processing  for  one-  and  two- 
dimensional  signals.  A  new  look  is  taken  at  multirate  system  theory  in  one  and  two 
dimensions  which  provides  the  framework  for  these  methodologies.  A  careful  analysis 
of  linear  optimal  filtering  for  problems  involving  different  input  and  output  sampling 
rates  is  conducted.  This  results  in  the  development  of  index  mapping  techniques  that 
simplify  the  formulation  of  Wiener-Hopf  equations  whose  solution  determine  the  op¬ 
timal  filters.  The  required  filters  exhibit  periodicity  in  both  one  and  two  dimensions, 
due  to  the  difference  in  sampling  rates.  The  reconstruction  algorithms  developed  are 
applied  to  one-  and  two-dimensional  reconstruction  problems. 
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EXECUTIVE  SUMMARY 


As  physical  and  manufacturing  limitations  are  reached  in  state-of-the-art  im¬ 
age  acquisition  systems,  there  is  increased  motivation  to  improve  the  resolution  of 
imagery  through  signal  processing  methods.  High-resolution  (HR)  imagery  is  desir¬ 
able  because  it  can  offer  more  detail  about  the  object  associated  with  the  imagery. 
The  “extra”  information  is  of  critical  importance  in  many  applications.  For  exam¬ 
ple,  HR  reconnaissance  images  can  provide  intelligence  analysts,  greater  information 
about  a  military  target,  including  its  capabilities,  operability  and  vulnerabilities,  and 
increase  analysts’  confidence  in  such  assessments.  Likewise,  HR  medical  images  can 
be  crucial  to  a  physician  in  making  a  proper  diagnosis  or  developing  a  suitable  treat¬ 
ment  regimen. 

Super-resolution  (SR)  image  reconstruction  is  an  approach  to  this  problem, 
and  this  area  of  research  encompasses  those  signal  processing  techniques  that  use 
multiple  low-resolution  (LR)  images  to  form  a  HR  image  of  some  related  object.  In 
this  work,  a  super-resolution  image  reconstruction  approach  is  proposed  from  the 
viewpoint  of  estimation  and  multirate  signal  processing  for  two-dimensional  signals 
or  images. 

Multirate  signal  processing  theory  deals  with  the  analysis  of  a  system  com¬ 
prised  of  multiple  signals  at  different  sampling  rates  and  is  fundamental  to  this  re¬ 
search.  An  example  of  such  a  system  is  a  sensor  network  that  collects  and  processes 
data  from  various  sensors,  where  the  information  from  each  sensor  might  be  collected 
at  a  different  rate.  In  developing  this  theory,  a  number  of  relationships  between  sig¬ 
nals  in  a  multirate  system  are  identified.  The  critical  finding  is  that  all  of  the  signals 
in  a  multirate  system  can  be  referred  to  a  single  “universal”  rate  for  that  system; 
therefore,  many  of  the  results  of  standard  signal  processing  theory  can  be  adapted  to 
multirate  systems  through  this  observation. 
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The  multirate  theory  developed  here  is  applied  to  signal  estimation,  where  one 
signal  is  estimated  from  some  other  related  signal  or  signals.  The  desired  signal  may 
be  corrupted  by  distortion  or  interference  and  is  usually  unobservable  (at  least  at 
the  moment  when  the  estimate  is  desired).  A  typical  signal  estimation  application  is 
the  recovery  of  a  transmitted  signal  from  a  received  signal  that  has  been  subject  to 
distortion  and  is  corrupted  by  noise. 

SR  image  reconstruction  can  be  viewed  as  a  problem  in  signal  estimation, 
where  a  related  LR  signal  or  signals  is  used  to  estimate  an  underlying  HR  signal. 
From  this  perspective,  the  observation  signal  or  signals,  and  desired  signal  form  a 
multirate  system.  This  motivates  the  application  of  the  theory  of  multirate  systems 
to  signal  estimation  and  the  resultant  extension  of  single-rate  signal  estimation  theory 
to  the  multirate  case. 

The  particular  branch  of  estimation  theory  applied  in  this  work  is  optimal 
filtering,  where  the  error  in  estimation  is  minimized  by  using  a  weighted  set  of  the 
LR  observation  images  to  filter  and  estimate  the  HR  image.  The  weights  used  in  this 
linear  estimate  are  called  filter  coefficients  and  application  of  this  theory  results  in 
a  set  of  equations  that  are  solved  to  obtain  these  coefficients  known  as  the  Wiener- 
Hopf  (WH)  equations.  In  this  research,  the  multirate  WH  equations  are  developed 
and  shown  to  have  a  periodically  time-dependent  solution.  Additionally,  the  concept 
of  index  mapping,  an  extension  of  the  multirate  theory,  is  developed  to  determine  the 
required  regions  of  the  LR  images  required  for  estimation. 

A  new  methodology  is  developed  and  presented,  by  application  and  extension 
of  the  results  of  multirate  and  optimal  estimation  theory  to  the  problem  of  SR  image 
reconstruction.  This  new  method  is  applied  to  a  set  of  LR  images,  and  the  resultant 
HR  image  is  compared  with  results  from  standard  interpolation  methods.  In  every 
case,  this  method  performed  better  than  the  standard  methods. 
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I.  INTRODUCTION 


As  physical  and  manufacturing  limitations  are  reached  in  state-of-the-art  im¬ 
age  acquisition  systems,  there  is  increased  motivation  to  improve  the  resolution  of 
imagery  through  signal  processing  methods.  Improvements  in  this  area  have  signifi¬ 
cant  commercial  and  military  application,  and  in  this  work  a  super-resolution  image 
reconstruction  approach  is  proposed  from  the  viewpoint  of  estimation  and  multirate 
signal  processing  for  two-dimensional  signals. 

A.  PROBLEM  STATEMENT/MOTIVATION 

Super-resolution  (SR)  imaging  has  recently  become  an  area  of  great  interest 
in  the  image  processing  research  community  (see  Section  I.B.2).  The  ability  to  form 
a  high-resolution  (HR)  image  from  a  collection  of  subsampled  images  has  a  broad 
range  of  applications  and  has  largely  been  motivated  by  physical  and  production 
limitations  on  existing  image  acquisition  systems  and  the  marginal  costs  associated 
with  increased  spatial  resolution.  Figure  1.1  depicts  the  SR  concept  where  a  collection 
of  low-resolution  (LR)  images  of  a  scene  are  superimposed  on  a  HR  grid,  available 
for  subsequent  HR  image  reconstruction. 

In  this  work,  we  propose  a  stochastic  multirate  approach  to  this  problem, 
adapting  and  extending  the  work  in  [Ref.  6,  7,  8,  9]  to  one-  and  two-dimensional 
signals.  The  earlier  work  has  focused  on  information  fusion  applications,  i.e.,  on  the 
combination  of  observations  from  multiple  sensors  to  perform  tracking,  surveillance, 
classification  or  some  other  task.  This  work  extends  these  concepts  to  reconstruction 
of  one-dimensional  signals  and  SR  image  reconstruction. 
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Figure  1.1.  Super-resolution  imaging  concept,  (After  [Ref.  1]). 


B.  PREVIOUS  WORK 

1.  Stochastic  Multirate  Signal  Processing 

Research  in  the  area  of  stochastic  multirate  signal  processing  has  been  lim¬ 
ited  to  a  handful  of  investigators  whose  work  has  focused  mainly  on  second  moment 
analysis  of  stochastic  systems,  from  both  temporal  and  spectral  points  of  view,  and 
optimal  estimation  theory,  including  both  Kalman  and  Weiner  filtering  theory. 

Vaidyanathan  et  al.  [Ref.  10,  11,  12]  investigate  how  the  statistical  properties 
of  stochastic  signals  are  altered  through  multirate  systems.  In  [Ref.  10] ,  several  facts 
and  theorems  are  presented  regarding  the  statistical  behavior  of  signals  as  they  are 
passed  through  decimators,  interpolators,  modulators,  and  more  complicated  inter- 
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connections.  For  example,  the  necessary  and  sufficient  condition  for  the  output  of 
an  L-fold  interpolation  filter  to  be  wide-sense  stationary  (WSS),  given  a  WSS  input, 
is  that  the  L-fold  decimation  of  the  filter  coefficients  results  in  no  aliasing,  i.  e.,  the 
filter  must  have  an  alias- free  (L)  support.  Additionally,  the  authors  illustrate  an 
application  of  this  theoretical  analysis  to  a  multirate  adaptive  filtering  scheme  for 
identification  of  band-limited  channels.  In  [Ref.  11],  this  work  is  continued  but  ad¬ 
dressed  using  bifrequency  maps  and  bispectra.  These  two-dimensional  (2-D)  Fourier 
transforms  characterize  all  linear  time- varying  (LTV)  systems  and  nonstationary  ran¬ 
dom  processes,  respectively.  In  fact,  by  using  these  concepts,  the  previous  results  are 
simplified  and  even  generalized  to  handle  the  case  of  vector  systems.  Finally,  in 
[Ref.  12],  further  analysis  is  conducted  using  bifrequency  maps  and  bispectra,  and  a 
bifrequency  characterization  of  lossless  LTV  systems  is  derived. 

Jahromi  et  al.  [Ref.  13,  14,  15]  consider  methods  to  optimally  estimate  samples 
of  a  random  signal  based  on  observations  made  by  multiple  observers  at  different 
sampling  rates  (lower  than  the  original  rate).  In  particular,  in  [Ref.  13],  the  problem 
of  fusing  two  low-rate  sensors  in  the  reconstruction  of  one  high-resolution  signal  is 
considered  when  time  delay  of  arrival  (TDOA)  is  present.  Using  the  “generalized 
cross-correlation”  technique,  the  delay  is  estimated  and  then  signal  reconstruction  is 
accomplished  using  perfect  reconstruction  synthesis  Liter  bank  theory.  In  [Ref.  14] 
and  [Ref.  15],  optimal  least  mean-square  estimation  is  used  to  develop  an  estimate 
for  samples  of  a  high-rate  signal.  The  estimator  is  a  function  of  the  power  spectral 
density  of  the  original  random  signal,  which  is  obtained  using  a  method  for  inductive 
inference  of  probability  distribution  referred  to  as  the  “maximum  entropy  principle” 
[Ref.  16], 

Chen  et  al.  [Ref.  17,  18,  19,  20]  investigate  use  of  the  Kalman  Liter  and 
Weiner  Liter  in  the  reconstruction  of  a  stochastic  signal  when  only  a  noisy,  downsam¬ 
pled  version  of  the  signal  can  be  measured.  In  [Ref.  17],  the  use  of  the  Kalman  Liter 
is  investigated  for  interpolating  and  estimating  values  of  an  autoregressive  or  moving 
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average  stochastic  signal  when  only  a  noisy,  downsampled  version  of  the  signal  can 
be  measured.  The  signal  reconstruction  problem  is  converted  into  a  state  estima¬ 
tion  problem  for  which  the  Kalman  filter  is  optimal.  Some  extensions  are  discussed, 
including  the  application  of  the  Kalman  reconstruction  filter  in  recovering  missing 
speech  packets  in  a  packet  switching  network  with  packet  interleaving.  Simulation 
results  are  presented,  which  indicate  that  the  multirate  Kalman  reconstruction  filters 
possess  better  reconstruction  performance  than  a  Wiener  reconstruction  filter  under 
comparable  numerical  complexity.  In  [Ref.  18] ,  a  multirate  deconvolution  filter  is  pro¬ 
posed  for  signal  reconstruction  in  multirate  systems  with  channel  noise.  Both  filter 
bank  and  transmultiplexer  architectures  are  used  to  demonstrate  the  design  proce¬ 
dure.  In  [Ref.  19],  a  block  state-space  model  is  introduced  where  transmultiplexer 
systems  unify  the  multirate  signals  and  channel  noise.  In  [Ref.  20],  the  optimal  signal 
reconstruction  problem  is  considered  in  transmultiplexer  systems  under  channel  noise 
from  the  viewpoint  of  Wiener-Hopf  theory.  A  calculus  of  variation  method  and  a 
spectral  factorization  technique  are  used  to  develop  an  appropriate  separation  filter 
bank  design. 

Scharf  et  al.  [Ref.  21]  introduce  a  least  squares  design  methodology  for  fil¬ 
tering  periodically  correlated  (PC)  scalar  time  series.  Since  any  PC  time  series  can 
be  represented  as  a  WSS  vector  time  series  where  each  constituent  subsequence  is 
a  decimated  version  of  the  original  shifted  in  time,  and  vice  versa,  multirate  filter 
banks  and  equivalent  polyphase  realizations  provide  a  natural  representation  for  this 
bidirectional  relationship.  This  relationship  affords  means  to  develop  a  spectral  rep¬ 
resentation  for  the  PC  time  series  and  hence  develop  causal  synthesis  and  causal 
whitening  filters  for  the  PC  scalar  time  series.  These  techniques  are  used  to  solve 
generalized  linear  minimum  mean-square  error  (MMSE)  filter  design  problems  for 
PC  scalar  time  series.  Note  that  this  viewpoint  can  be  extended  to  multirate  systems 
where  the  correlation  between  observation  sequences  is  periodically  correlated. 
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Therrien  et  al.  [Ref.  6,  22,  7,  8,  9,  23]  develop  theory  and  methodology 
required  for  employing  optimal  linear  filtering  in  estimating  an  underlying  signal 
from  observation  sequences  at  different  sampling  rates.  The  focus  of  these  efforts  is 
on  information  fusion,  i.e. ,  on  the  combination  of  observations  from  multiple  sensors 
to  perform  tracking,  surveillance,  classification  or  some  other  task.  In  particular, 
[Ref.  6],  [Ref.  22]  and  [Ref.  7]  consider  a  simplified  problem  where  an  underlying 
signal  is  estimated  from  two  sequences,  one  observed  at  full  rate  and  the  other  at 
half  the  rate.  In  [Ref.  8],  least  squares  formulations  are  examined  where  the  second 
sequence  has  an  arbitrary  sampling  rate.  In  [Ref.  9],  a  general  approach  is  suggested 
for  any  number  of  observation  signals  at  arbitrary  sampling  rates.  Finally,  in  [Ref. 
23],  previous  theory  and  methods  are  developed  to  consider  the  problem  of  HR  signal 
and  image  reconstruction.  This  work  forms  the  basis  for  the  proposed  research  and 
represents  an  advance  in  the  area  of  super-resolution  image  reconstruction. 

2.  Super-Resolution  Reconstruction/Imaging 

Generally,  super-resolution  (SR)  image  reconstruction  refers  to  signal  process¬ 
ing  methods  in  which  a  high-resolution  (HR)  image  is  obtained  from  a  set  or  ensemble 
of  observed  low-resolution  (LR)  images  [Ref.  1],  If  each  observed  LR  image  is  sub¬ 
sampled  (and  aliased)  and  is  translated  by  a  different  subpixcl  amount,  this  set  of 
unique  observation  images  can  be  used  for  reconstruction.  Figure  1.1  demonstrates 
this  conceptually.  Both  [Ref.  1]  and  [Ref.  2]  provide  general  surveys  of  research  to 
date  regarding  this  topic,  and  the  following  major  areas  of  research  are  identified: 
nonuniform  interpolation,  frequency  domain,  regularized  SR  reconstruction,  projec¬ 
tion  onto  convex  sets  (POCS),  maximum  likelihood  (ML)  projection  onto  convex  sets 
(ML-POCS)  hybrid  reconstruction,  and  other  approaches  [Ref.  1], 

The  most  prevalent  approaches  in  the  literature  are  those  based  on  nonuni¬ 
form  interpolation.  These  approaches  typically  use  a  three-stage  sequential  process, 
comprised  of  registration,  interpolation,  and  restoration.  The  registration  step  is  a 
mapping  of  pixels  from  each  LR  image  to  a  reference  grid,  which  results  in  a  HR  grid 
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comprised  of  a  set  of  nonuniformly  spaced  pixels.  The  interpolation  step  conforms 
these  nonuniformly  spaced  pixels  to  a  uniform  sampling  grid,  which  results  in  the 
upsampled  HR  image.  Finally,  restoration  removes  the  effects  of  sensor  distortion 
and  noise.  This  scheme  is  depicted  in  Figure  1.2.  Representative  works  include  [Ref. 
24,  25,  26,  27], 


LR  Images 


HR  Image 

y 


Figure  1.2.  Typical  model  for  nonuniform  interpolation  approach  to  SR,  (From  [Ref. 

2])- 


The  frequency-domain  approaches  exploit  the  relationship  between  the  discrete 
Fourier  transforms  (DFT)  of  the  LR  images  and  the  continuous  Fourier  transform 
(CFT)  of  the  desired  HR  image  by  using  the  information  generated  through  relative 
motion  between  the  LR  images,  the  aliasing  generated  by  downsampling  relative  to 
the  desired  HR  image,  and  the  assumption  that  the  original  HR  image  is  bandlim- 
ited.  A  set  of  linear  system  equations  are  developed,  and  the  continuous  Fourier 
coefficients  are  found.  The  desired  HR  image  is  estimated  from  the  CFT  synthesis 
equation.  Tsai  and  Huang  [Ref.  28]  were  the  first  to  introduce  this  method  and  were 
also  the  first  researchers  to  address  the  problem  of  reconstructing  a  HR  image  from  a 
set  of  translated  LR  images.  Kim  et  al.  [Ref.  29]  extended  this  approach  to  include 
the  presence  of  noise  in  the  LR  images  using  a  recursive  procedure  based  on  weighted 
least  squares  theory.  Kim  and  Su  [Ref.  30]  further  extended  this  approach  by  consid- 
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ering  noise  and  different  blur  distortions  in  the  LR  images.  Vandewalle  et  al.  [Ref. 
31]  consider  offset  estimation  using  a  subspace  minimization  method  followed  by  a 
frequency-based  reconstruction  method  based  on  the  continuous  and  discrete  Fourier 
series. 

The  regularized  SR  reconstruction  methods  use  regularization  methods  to  solve 
the  often  ill-posed  inverse  problem  introduced  in  the  frequency-domain  approaches. 
Typically,  the  ill-posed  problems  are  a  result  of  an  insufficient  number  of  LR  images 
or  ill-conditioned  blur  operators  [Ref.  1],  Generally,  two  approaches  have  been  con¬ 
sidered:  deterministic  and  stochastic  regularization.  Deterministic  approaches  [Ref. 
32,  33,  34,  35]  typically  use  constrained  least  squares  methods  (CLS)  while  stochastic 
approaches  [Ref.  36,  37,  38]  typically  use  maximum  a  posteriori  (MAP)  or  maximum 
likelihood  (ML)  methods. 

POCS  methods  are  based  on  set  theoretic  estimation  theory  [Ref.  39].  Rather 
than  using  conventional  estimation  theory,  the  POCS  formulations  incorporate  a  pri¬ 
ori  knowledge  into  the  solution  and  yield  a  solution  consistent  with  user-furnished 
constraints.  Application  of  this  method  as  applied  to  SR  was  introduced  by  Stark 
and  Oskoui  [Ref.  40]  and  extended  by  Tekalp  et  al.  in  [Ref.  41,  42],  which  takes 
into  account  the  presence  of  both  sensor  blurring  and  observation  noise,  and  suggests 
POCS  as  a  new  method  for  restoration  of  spatially-variant  blurred  images. 

ML-POCS  hybrid  reconstruction  approaches  estimate  desired  HR  images  by 
minimizing  the  ML  or  MAP  cost  functional  while  constraining  the  solution  within 
certain  closed  convex  sets  in  accordance  with  POCS  methodology  [Ref.  37]. 

There  are  a  number  of  other  areas  that  are  considered  in  the  literature,  and 
some  examples  are  presented  here.  One  approach  attempts  to  reconstruct  a  HR  image 
from  a  single  LR  image  and  is  referred  to  as  improved  definition  image  interpolation 
[Ref.  43].  Another  area  of  study,  referred  to  as  iterative  back-projection  [Ref.  44,  45, 
46],  uses  tomographic  projection  methods  to  estimate  a  HR  image.  Researchers  are 
also  considering  the  SR  problem  when  no  relative  subpixel  motion  exists  between  LR 
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images.  By  considering  differently  blurred  LR  images,  motionless  SR  reconstruction 
can  be  demonstrated  [Ref.  47,  48].  Milanfar  et  al.  analyze  the  joint  problem  of 
image  registration  and  HR  reconstruction  in  the  context  of  fundamental  statistical 
performance  limits.  By  using  the  Cramer-Rao  bound,  they  demonstrate  ability  to 
bound  estimator  performance  in  terms  of  MSE,  examining  performance  limits  as 
they  relate  to  such  imaging  system  parameters  as  the  downsampling  factor,  signal-to- 
noise  ratio,  and  point  spread  function.  Finally,  researchers  are  considering  adaptive 
filtering  approaches  to  the  SR  problem,  considering  modified  recursive  least  squares 
(RLS),  linear  mean-square  (LMS)  and  steepest  descent  methods  [Ref.  49]. 

C.  THESIS  ORGANIZATION 

This  manuscript  is  organized  as  follows.  The  current  chapter  is  introductory 
and  presents  the  motivation  for  this  work,  defining  the  problem  and  outlining  the 
approach  used  to  solve  it.  Additionally,  a  review  of  the  relevant  literature  is  included, 
both  in  the  area  of  stochastic  multirate  signal  processing  and  super-resolution  image 
reconstruction. 

The  second  chapter  introduces  various  fundamental  signal  processing  and 
mathematical  concepts  required  for  theoretic  and  application-related  developments 
in  future  chapters.  These  include  various  signal  taxonomies  and  representations,  a 
review  of  relevant  topics  in  second-moment  analysis,  and  required  number  theory  and 
linear  algebra  concepts.  Further,  this  chapter,  establishes  notation  and  conventions 
for  purposes  of  consistency  throughout  this  work. 

In  the  third  chapter,  the  theory  of  multirate  systems  is  established.  In  this 
analysis,  the  relationships  between  a  multirate  system  and  its  constituent  signals  are 
characterized,  the  system  theory  for  multirate  systems  is  developed,  and  the 


representation  of  discrete  linear  systems  is  presented  from  a  system  theoretic  point 
of  view.  Finally,  a  linear  algebraic  approach  is  introduced  to  model  various  multirate 
operations  for  use  in  reconstruction  applications. 

Chapter  IV  develops  the  concept  of  multirate  signal  estimation  and  is  founda¬ 
tional  in  developing  stochastic  approaches  to  solving  the  signal  reconstruction  prob¬ 
lem.  The  optimal  filtering  problem  is  introduced  in  terms  of  the  ordinary  Wiener- 
Hopf  equation  and  is  then  expanded,  first  to  the  single-channel,  multirate  estimation 
problem  and  then  to  the  multi-channel,  multirate  problem.  Also  in  this  chapter,  the 
relationship  between  samples  in  one  signal  domain  to  those  in  a  different  signal  do¬ 
main  (signals  at  different  rates)  is  established  through  the  concept  of  index  mapping, 
which  allows  for  a  very  general  representation  of  the  multirate  Wiener-Hopf  equations. 

Chapter  V  considers  the  problem  of  signal  reconstruction  in  the  one-  and  two- 
dimensions.  In  this  chapter,  the  problem  is  stated  for  both  cases,  observation  models 
are  established,  reconstruction  approaches  and  algorithms  are  developed,  and  then 
the  results  of  each  algorithm  are  presented. 

Finally,  Chapter  VI  provides  conclusory  remarks  on  the  hirelings  of  this  re¬ 
search  and  establishes  direction  for  future  work  related  to  this  research. 
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II.  PRELIMINARIES,  CONVENTIONS,  AND 

NOTATION 

In  the  development  of  approaches  to  signal  and  image  reconstruction,  a  num¬ 
ber  of  fundamental  concepts  from  the  areas  of  signal  processing  and  mathematics  are 
required.  In  this  chapter,  a  foundation  is  set  in  these  areas  upon  which  the  theory  of 
multirate  signals  and  multirate  estimation  will  be  built.  In  doing  so,  we  present  the 
underlying  concepts,  but  also  emphasize  required  definitions,  notations  and  conven¬ 
tions,  in  order  to  ensure  consistency  and  accuracy,  and  to  facilitate  understanding. 

A.  SIGNALS 

1.  Etymology 

Etymologically  speaking,  the  word  signal  is  derived  from  the  Latin  signum, 
which  can  be  rendered  as  “a  sign,  mark,  or  token;”  or  in  a  military  sense,  “a  standard, 
banner,  or  ensign;”  or  “a  physical  representation  of  a  person  or  thing,  like  a  figure, 
image,  or  statue  [Ref.  50].”  Generally,  the  Latin  seems  to  imply  that  a  signum  is 
something  that  conveys  information  about  or  from  someone  or  something  else.  The 
relevant  modern  dictionary  definition  of  signal  carries  this  idea  further:  “a  detectable 
physical  quantity  or  impulse  by  which  messages  or  information  can  be  transmitted 
[Ref.  51].” 

In  the  area  of  electrical  engineering  known  as  digital  signal  processing,  a  related 
but  more  helpful  definition  of  a  signal  is  a  collection  of  information,  usually  a  pattern 
of  variation  [Ref.  52],  that  describes  some  physical  phenomenon.  In  other  words,  a 
signal  conveys  relevant  information  about  some  physical  phenomena  (signum).  The 
variation  in  electrical  voltage  measured  at  the  input  of  an  electronic  circuit,  the 
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variation  in  acoustic  pressure  sensed  by  a  microphone  recording  a  musical  concert, 
or  the  variation  in  light  intensity  captured  by  a  camera  recording  a  scene  are  all 
examples  of  signals  treated  in  modern  signal  processing. 

2.  Signal  Definitions 

Throughout  this  presentation,  various  types  of  signals  and  sequences  are  in¬ 
troduced  and  analyzed.  In  this  section,  for  the  sake  of  clarity,  the  definition  of  such 
signals  and  sequences  are  established,  as  are  the  associated  conventions  and  nota¬ 
tions.  Let  us  begin  with  one-dimensional  signals  that  are  scalar-valued.  We  define 
these  more  precisely  below. 

a.  Deterministic  Signals  and  Sequences 
A  deterministic  analog  signal  or  simply  an  analog  signal  is  defined  as 

follows. 

Definition  1.  A  deterministic  analog  signal ,  denoted  by  (a;(t)},  or  when  it  is  clear 
from  context  x(t),  is  a  set  of  ordered  measurements  such  that  for  every  t  G  M,  there 
exists  a  corresponding  measurement  m  =  x(t).  If  all  such  measurements  are  members 
of  the  extended  real  numbers1,  then  x{t)  is  said  to  be  a  real-valued  (or  real )  analog 
signal.  If  the  measurements  are  members  of  the  complex  numbers,  then  the  signal  is 
said  to  be  a  complex-valued  (or  complex )  analog  signal. 

An  analog  signal  is  frequently  represented  by  a  mathematical  function, 
which  may  or  may  not  be  continuous.  For  example,  the  signal  known  as  the  unit-step , 
defined  by 

{1  t  >  0 

(2.1) 

0  t  <  0 

is  well  known  in  signal  processing,  but  the  function  representing  it  is  not  continuous 
(at  t  —  0). 

1The  extended  real  numbers  are  defined  as  R  =  1  U  {- oo,oo}. 
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Although  many  signals  are  represented  by  functions  defined  on  the  real 
number  line,  our  definition  of  a  signal  is  not  necessarily  the  same  as  the  mathematical 
definition  of  a  function.  The  set  of  analog  signals  commonly  includes  the  unit  impulse, 
which  (strictly  speaking)  is  not  a  function  at  all  but  a  distribution  or  “generalized 
function,”  described  by  a  careful  limiting  process  [Ref.  53,  54]  to  insure  that  the 
resulting  entity  satisfies  certain  conditions  when  it  appears  in  an  integral. 

Signals  may  have  many  other  properties  that  provide  for  further  char¬ 
acterization.  One  property  of  concern  in  this  work  is  that  of  periodicity.  A  signal  is 
said  to  be  periodic  if  there  exists  a  positive  real  number  P  such  that 

x(t )  =  x(t  +  P )  for  all  t.  (2.2) 

The  smallest  such  P  is  called  the  period. 

A  deterministic  sequence  (or  simply  a  sequence)  is  defined  as  follows. 

Definition  2.  A  deterministic  sequence,  denoted  by  {x[n]},  or  when  clear  from  con¬ 
text  x [n] ,  is  a  countable  set  of  ordered  measurements  such  that  for  every  n  e  Z,  there 
exists  a  corresponding  measurement  m  —  x  [n] .  If  all  such  measurements  are  members 
of  the  extended  real  numbers,  then  x [n  is  said  to  be  a  real-valued  (or  real)  sequence. 
If  the  measurements  are  members  of  the  complex  numbers,  then  the  sequence  is  said 
to  be  a  complex-valued  (or  complex)  sequence. 

A  sequence  x[n]  is  said  to  be  periodic  if  there  exists  a  positive  integer 

N  such  that 

x[n]  =  x[n  +  N]  for  all  n,  (2.3) 

and  the  smallest  such  N  is  called  the  period.  Note  that  not  all  sequences  derived 
by  sampling  a  periodic  analog  signal  are  periodic.  For  example,  the  analog  signal 
x(t)  =  cos(27t fot  +  (f>)  is  periodic  for  any  real  number  /0,  while  the  sequence  x[n] 
defined  by  x[n]  =  x(nTs)  =  cos(2nfonTs  +  0)  is  periodic  only  if  f0Ts  is  a  rational 
number. 
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Observe  that  both  a  signal  and  a  sequence  are  defined  by  an  ordered 
set  of  measurements,  but  over  a  different  domain  (R  or  Z).  Further,  parentheses  are 
used  in  the  notation  for  an  analog  signal  x(-)  while  square  brackets  are  used  for  a 
sequence  &[•]  (to  indicate  the  discrete  nature  of  its  domain).  The  variable  t  or  n  is 
frequently  used  to  represent  time,  although  the  units  of  “time”  need  to  be  specified 
in  any  real-world  problem.  In  the  case  of  a  sequence,  n  is  just  an  index  variable  used 
to  order  the  measurements,  and  there  is  need  in  signal  processing  to  define  what  will 
be  called  a  deterministic  discrete- domain  signal  or  simply  discrete-domain  signal  . 

Definition  3.  A  deterministic  discrete- domain  signal,  denoted  by  {xxit)},  or  when 
it  is  clear  from  context  Xt (t),  is  a  set  of  ordered  measurements  such  that  for  every 
f  Efr,  there  exists  a  corresponding  measurement  m  =  xr(t),  where 
d'r  =  {' nT ;  n  G  Z},  and  T  is  a  positive  real  number  called  the  sampling  interval.  The 
signal  domain  is  defined  as  the  set  x1jt-  If  all  such  measurements  are  members  of  the 
extended  real  numbers,  then  xr(t)  is  said  to  be  a  real-valued  (or  real)  discrete-domain 
signal.  If  the  measurements  are  members  of  the  complex  numbers,  then  the  signal  is 
said  to  be  a  complex-valued  (or  complex)  discrete-domain  signal.  When  t  represents 
time,  a  discrete-domain  signal  may  be  called  a  discrete-time  signal. 

This  definition  of  a  discrete-domain  signal  is  similar  to  that  of  an  analog 
signal  except  that  the  signal  is  defined  on  a  countable  set  An  important  obser¬ 
vation  is  that  a  discrete-domain  signal  is  equivalent  to  a  sequence  and  an  associated 
sampling  interval  T  or  its  reciprocal  F  —  1/T, 

xx(t)  =  {x[n\,T}  =  {x[n],F}  for  n  e  Z.  (2.4) 

The  quantity  F  is  called  the  sampling  rate  (in  samples/sec  or  Hz)  and  in  discussing 
discrete-domain  signals,  it  is  common  to  refer  to  the  sequence  and  its  sampling  rate. 
For  example,  we  may  use  the  expression  “x [n]  at  a  rate  of  20  kHz”  to  describe  a 
discrete-domain  signal,  which  has  a  sampling  interval  of  T  =  0.05  msec. 
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It  is  also  common  not  to  mention  the  sampling  rate  if  the  sampling 
rate  is  common  throughout  a  system  (single-rate  system).  On  the  other  hand,  when 
dealing  with  a  multirate  system,  it  is  common  to  use  different  letters,  such  as  n  and  m, 
to  designate  sequences,  for  example,  x[n)  and  y [m] ,  to  indicate  that  these  sequences 
represent  discrete-domain  signals  with  different  sampling  rates. 

Figure  2.1  illustrates  a  discrete-domain  signal.  Note  that  the  signal  is 
defined  only  on  the  points  t  =  nT  and  is  undefined  everywhere  else.  Note,  also,  that 
while  a  discrete-domain  signal  may  be  derived  by  sampling  an  analog  signal,  this  is  not 
always  the  case.  Any  sequence,  regardless  of  how  it  is  computed  (say  in  MATLAB  or 
on  an  ASIC  chip)  when  combined  with  a  sampling  interval,  defines  a  discrete-domain 
signal.  The  corresponding  analog  signal  need  not  exist  unless  (as  in  the  output  of  a 
digital  signal  processing  chain)  some  special  action  is  taken  to  construct  it. 


xT(t) 


Figure  2.1.  Graphical  representation  of  a  discrete-domain  signal  Xr(t)  with  sampling 
interval  T  =  0.05.  Note  that  the  signal  is  defined  only  at  t  =  nT ;  n  6  Z. 

b.  Random  Signals  and  Sequences 

In  statistical  signal  processing,  a  probabilistic  model  is  necessary  for 
signals.  This  model  is  embedded  in  the  concept  of  a  random  signal  or  a  stochastic 
signal.  A  real  random  signal  or  ( real  stochastic  signal )  is  defined  as  follows. 
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Definition  4.  A  real  random  signal ,  denoted  by  (A"(f)},  or  when  it  is  clear  from 
context  X(t),  is  a  set  of  ordered  random  variables  (representing  measurements)  such 
that  for  every  tel,  there  exists  a  corresponding  random  variable  X(t). 

Note  that  when  the  context  is  clear,  a  random  signal  may  be  designated  by  a  lower 
case  variable,  i.e.,  x(t),  d(t),  etc. 

Since  a  random  variable  is  a  mapping  from  some  sample  space  to  the 
real  line,  the  definition  for  a  complex  random  signal  requires  special  caution.  The 
following  definition  is  therefore  provided. 

Definition  5.  A  complex  random  signal  or  ( complex  stochastic  signal ),  denoted  by 
(Z(t)},  is  defined  by  Z(t)  =  X(t)+jY(t),  where  X(t)  and  Y(t )  are  real  random  analog 
signals  defined  on  a  common  domain.  In  other  words,  for  every  tel,  there  exists  a 
pair  of  corresponding  random  variables  X(t)  and  Y(t)  such  that  Z(t)  =  X (t)  +jY(t). 
Again,  we  may  use  Z(t)  instead  of  {Z(t)}  when  the  meaning  is  clear  from  context. 

Random  sequences  and  random  discrete-domain  signals  can  be  defined 
in  a  similar  manner. 

Definition  6.  A  real  random  sequence  or  ( real  stochastic  sequence ),  denoted  by 
{A[n]},  is  a  countable  set  of  ordered  random  variables  (representing  measurements) 
such  that  for  every  n  G  Z,  there  exists  a  corresponding  random  variable  X[n].  A 
complex  random  sequence  can  be  defined  in  a  manner  similar  to  that  of  a  complex 
random  signal. 

Note  that  when  the  context  is  clear,  a  random  sequence  may  be  designated  by  a  lower 
case  variable,  i.e.,  x[n],  d[n\,  etc. 

Definition  7.  A  random  discrete-domain  signal,  denoted  by  {A A(f)},  or  when  it  is 
clear  from  context  Xr(t),  is  a  set  of  ordered  random  variables  (representing  mea¬ 
surements)  such  that  for  every  t  G  xI>t,  there  exists  a  corresponding  random  variable 
Afr(t),  where  =  {nT;n  G  Z},  and  T  is  the  sampling  interval. 


16 


A  random  discrete-domain  signal  is  sometimes  also  referred  to  as  a  time  series ;  how¬ 
ever,  the  use  of  that  term  in  the  literature  is  not  always  consistent. 

c.  Multi-channel  Signals  and  Sequences 
In  signal  processing,  it  is  often  the  case  that  a  system  may  contain 
signals  or  sequences  that  are  derived  from  multiple  sources  or  multiple  sensors.  In 
order  to  represent  such  signals  and  sequences,  multi-channel  signals  and  sequences 
are  defined.  A  multi-channel  signal  is  a  set  of  (single-channel)  signals  that  share  a 
common  domain  and  is  represented  by  a  vector 

aq(f)"| 

^2 

XN(t)  J 

whose  components  Xi(t),  x2(t), . . . ,  apv(t)  are  (analog  or  discrete-domain)  signals  as 
defined  earlier.  The  signals  may  be  real  or  complex,  deterministic  or  random.  By 
convention,  bold  face  and  vector  notation  are  used  to  represent  such  signals  as  in 


x(f)  = 


x(f) 


cosc at 
—  sin  cut 


or  m 


X(i)  = 


A  cos(ut  +  <E>) 

—A  sin(c at  +  $) 

where  X(t )  represents  a  random  signal  defined  by  random  variables  A  and  $. 
A  multi-channel  sequence 


x[n 


x\  [n\ 
x2[n\ 


Xn[u 
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is  represented  by  a  vector  whose  components  x i  [n] ,  x2 [n] , . . .  ,xn[ti\  are  sequences  as 
defined  earlier.  Again,  all  of  the  terms  describing  an  individual  sequence  (e.g.,  real, 
complex,  etc.)  can  be  applied  to  a  multi-channel  sequence. 

d.  Two-dimensional  Signals  and  Sequences 

Since  two-dimensional  signals  and  sequences  are  at  the  heart  of  image 
processing,  it  is  helpful  to  characterize  the  2-D  counterparts  to  the  familiar  one¬ 
dimensional  signals  and  sequences  already  presented.  A  two-dimensional  (2-D)  analog 
signal  is  defined  as  follows. 

Definition  8.  A  two-dimensional  (2-D)  analog  signal,  denoted  by  {x(t\,  £2)},  or  when 
it  is  clear  from  context  x(ti,t2),  is  a  set  of  ordered  measurements  such  that  for  every 
pair  (f  1, f2)  €  M2,  there  exists  a  corresponding  measurement  m  =  x(fi,f2)-  Two- 
dimensional  signals  can  be  real  or  complex,  deterministic  or  random,  ft  is  sometimes 
convenient  to  represent  a  2-D  signal  with  a  bold  face  argument  t  =  (ti,  f2)  G  M2.  Thus, 
the  2-D  signal  would  be  denoted  by  {x(t)}  or  x(t)  when  clear  from  the  context. 

Although  a  sequence  seems  to  imply  an  ordered  set  of  terms  in  one 
dimension,  it  is  common  in  signal  processing  to  extend  the  meaning  to  apply  to 
signal  defined  on  a  two-dimensional  domain.  A  two-dimensional  sequence  and  two- 
dimensional  discrete-domain  signal  are  thus  defined  as  follows. 

Definition  9.  A  two-dimensional  sequence,  denoted  by  {x[ni,n2]},  or  when  it  is 
clear  from  context  x[n\,  n2],  is  a  set  of  ordered  measurements  such  that  for  every  pair 
(711,  n2)  G  Z2,  there  exists  a  corresponding  measurement  m  =  x(ri\,  n2].  2-D  sequences 
can  be  real  or  complex,  deterministic  or  random;  they  may  also  be  represented  as 
{a:[n]}  or  x[n],  where  the  boldface  argument  denotes  the  ordered  pair  {ri\,  n2)  e  Z2. 

Definition  10.  A  two-dimensional  discrete- domain  signal,  denoted  by  {^r1r2(D,  h)} 
or  xr1r2(f  1,  t2),  is  a  set  of  ordered  measurements  such  that  for  every  pair  {t\,  f2)  in  the 
domain  vI;t1t2  =  x  where  is  as  defined  earlier,  there  exists  a  corresponding 
measurement  m  =  xt1t2(D,  t2),  and  Ti  and  T2  are  the  associated  sampling  intervals. 
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For  convenience  in  notation,  we  may  use  XT(t)  and  'Ft  to  denote  the  2-D  signal 
and  its  domain,  where  T  represents  the  ordered  pair  (Ti,!^)  of  sampling  intervals. 
Again,  note  that  a  two-dimensional  discrete-domain  signal  can  be  real  or  complex, 
deterministic  or  random. 

The  image  projected  on  the  him  plane  of  a  camera  is  an  example  of 
a  2-D  analog  signal.  If  him  is  thought  of  as  a  continuous  medium,  then  the  image 
captured  on  the  him  is  also  a  representation  of  a  2-D  analog  signal.  If  the  image  is 
projected  onto  a  sensor  array  as  in  a  digital  camera,  then  the  resulting  sampled  image 
is  represented  by  a  2-D  discrete-domain  signal. 

Signals  can  be  both  multi-dimensional  and  multi-channel.  A  common 
example  is  a  color  image  where  the  domain  is  two-dimensional  (horizontal  and  vertical 
spatial  variables),  and  there  are  3  channels  corresponding  to  the  three  components  of 
a  color  space,  such  as  RGB  (red,  green,  blue),  CMY  (cyan,  magenta,  yellow)  or  HSI 
(hue,  saturation,  intensity). 

Two-dimensional  random  signals  and  sequences  are  similar  to  their  cor¬ 
responding  deterministic  representations  except  that  the  measurements  are  repre¬ 
sented  by  random  variables. 
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e.  Summary  of  Notation  and  Convention 

A  summary  of  the  various  signal  representations  is  provided  in  Ta¬ 
ble  2.1. 


Representation 

Name 

x{t) 

Deterministic  analog  signal,  analog  aignal 

x[n] 

Deterministic  sequence 

xT(t),x[n]T 

Deterministic  discrete-domain  signal  with  sampling  in¬ 
terval  T,  Discrete-domain  signal 

x(t1,t2),x(t) 

Two-dimensional  deterministic  analog  signal,  2-D  analog 
signal 

x[ni,n2],x[n] 

Two-dimensional  deterministic  sequence,  2-D  determin¬ 
istic  sequence 

XTlT2(tl,  t2),  2>r(t) 

Two-dimensional  deterministic  discrete-domain  signal 
with  sampling  intervals  T\  and  T2l  2-D  discrete-domain 
signal 

X{t) 

Random  analog  signal 

X[n\ 

Random  sequence 

Table  2.1.  Summary  of  signal  representations. 


B.  CONCEPTS  IN  LINEAR  ALGEBRA 

1.  Random  Vectors 

Often,  it  is  necessary  to  process  some  finite  number  of  samples  of  a  random 
sequence.  Such  a  finite-length  sequence  can  be  conveniently  represented  by  a  random 
vector  [Ref.  5].  This  provides  for  compact  notation  and  formulation  and  solution  of 
problems  in  a  linear  algebra  sense.  A  random  sequence  X[n]  restricted  to  some 
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interval  0  <  n  <  N  —  1  can  be  represented  by  an  TV-component  random  vector  x  as 
shown  in  Figure  2.2  and  written  as 


.Y[0] 


Figure  2.2.  Graphical  representation  of  a  finite-length  random  sequence  as  a  random 
vector. 


2.  Kronecker  Products 

The  Kronecker  product,  also  known  as  the  direct  product  or  tensor  product, 
has  its  origins  in  group  theory  [Ref.  4]  and  has  important  applications  in  a  number  of 
technical  disciplines.  In  this  study,  the  Kronecker  product  is  used  to  develop  matrix 
representations  of  various  multirate  operations. 


Definition  11.  Let  A  be  an  m  x  n  matrix  (with  entries  a^)  and  let  B  be  an  r  x  s 
matrix.  Then  the  Kronecker  product  of  A  and  B  is  the  mr  x  ns  block  matrix 

i\ 


A®  B  = 


f  onB 

ai2B  . 

a2iB 

a22B 

•  ^2  n  B 

yamiB 

®m2B 

(2,6) 


Equation  (2.6)  is  also  called  a  right  Kronecker  product  as  opposed  to  the 
definition  A  (g/  B  =  B  ®  A,  which  is  called  a  left  Kronecker  product.  Since  there  is 
no  need  to  use  both,  we  will  stick  with  the  more  common  definition  (2.6). 
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A  summary  of  some  important  properties  of  the  Kronecker  product  is  provided 


in  Table  2.2. 


A  (g)  (aB)  =  a(A  <8>  B) 

(A  +  B)<g>C  =  A<g)C  +  B<8)C 
A  ®  (B  ®  C)  =  (A  (8)  B)  (8)  C 
(A®B)t  =  At®Bt 
(A  <g>  B)(C  <g>  D)  =  AC  <g>  BD 
(A  OB)’1  =  A'1  OB-1 


Table  2.2.  Some  Kronecker  product  properties  and  rules,  (After  [Ref.  4]). 


3.  Reversal  of  Matrices  and  Vectors 

In  signal  processing,  it  is  a  common  requirement  to  view  signals  as  evolving 
either  forward  or  backward  in  time.  A  well-known  example  is  the  convolution  opera¬ 
tion,  where  the  linear  combination  of  terms  involves  a  time-reversed  version  of  either 
the  input  signal  or  the  system  impulse  response.  Since,  in  discrete-time  signal  pro¬ 
cessing,  signals  are  often  represented  by  vectors,  it  is  useful  to  define  the  operation 
of  reversal  for  vectors  and  matrices. 

The  reversal  of  a  vector  x  is  the  vector  with  its  elements  in  reverse  order. 
Given  the  vector 


Xi 

XN 

X  = 

X2 

,  its  reversal  is  x  = 

Xn-1 

XN 

Xi 

Note  that  the  notation  for  the  reversal  is  x,  and  it  is  used  just  like  notation  for  the 
transposition  of  a  vector  or  matrix. 


22 


The  reversal  of  a  matrix  A  is  the  matrix  with  its  column  and  row  elements  in 


reverse  order.  Given  the  matrix  A  G  M. 

MxN 

an 

ai2 

OlN 

A  = 

a2i 

0-22 

02N 

Ami 

«M  2 

omn 

its  reversal  A  G  MMxiV  is  given  by 

OMN 

«M2 

Om  1 

A  = 

O2N 

0-22 

021 

a±N 

Ol2 

Oil 

Note  that  the  reversal  of  a  vector  or  matrix  can  be  formed  by  the  product  of  a 
conformable  counter  identity  and  the  vector  or  matrix  itself. 

Some  common  properties  of  the  reversal  operator  are  included  in  Table  2.3. 
In  particular,  the  reversal  of  matrix  and  Kronecker  products  (see  Section  II. B. 2)  are 
products  of  the  reversals,  and  the  operation  of  reversal  commutes  with  inversion, 
conjugation  and  transposition. 


Quantity 

Reversal 

Matrix  product 

AB 

AB 

Matrix  inverse 

A-1 

(A)’1 

Matrix  conjugate 

A* 

(A)* 

Matrix  transpose 

AT 

(Af 

Kronecker  product 

A  ®  B 

A  ®  B 

Table  2.3.  Some  properties  of  the  reversal  operator,  (After  [Ref.  5]). 
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4.  Frobenius  Inner  Product 

In  the  development  of  approaches  to  two-dimensional  signal  reconstruction,  it 
is  convenient  to  express  the  related  linear  estimates  in  terms  of  the  Frobenius  inner 
product. 

Definition  12.  For  any  A,  B  e  Wnxn,  with  elements  a%3 ,  bl3 ,  the  Frobenius  inner 
product  of  the  matrices  is  defined  as 

m  n 

(A,  B)  =  tr(ABr)  =  ££«*<■  <*») 

i=i  j= i 

C.  MOMENT  ANALYSIS  OF  RANDOM  PROCESSES 

Generally,  a  complete  statistical  model  is  unavailable  when  analyzing  systems 
of  random  processes.  Either  the  required  joint  density  functions  are  unavailable,  or 
they  are  too  complex  to  be  of  utility.  If  the  random  processes  under  consideration  are 
Gaussian,  then  the  system  can  be  fully  specified  by  only  its  first  two  moments  [Ref.  5] . 
Even  if  the  processes  are  not  Gaussian,  second  moment  analysis  is  often  adequate  in 
characterizing  the  statistical  relationships  between  signals  in  such  systems  and  forms 
the  basis  for  any  additional  analyses.  This  section  introduces  the  required  definitions 
and  relevant  properties  associated  with  second  moment  analysis  [Ref.  5]. 

1.  Definitions  and  Properties 

Given  the  random  process  X [n] ,  the  first  moment  or  mean  of  the  random 
process  is  defined  by 

mx[n]  =  £{X[n\},  (2.10) 

where  £{•}  denotes  expectation. 

The  correlation  between  any  two  samples  of  the  random  process  X[ni]  and 
A  [no]  is  described  by  the  correlation  function  or  autocorrelation  function,  which  is 
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defined  by 


Rx  [ni,n0\  =  £{X[ni]X*  [n0]}.  (2.11) 

In  certain  applications,  and  extensively  in  this  work,  it  is  convenient  to  define 
a  time- dependent  correlation  function  as 

Rx[n\  l }  =  E{X[n]X*[n  -  /]},  (2.12) 

and  the  various  definitions  and  relationships  introduced  in  this  section  will  be  based 
on  this  “time-dependent”  representation. 

The  covariance  between  any  two  samples  of  the  random  process  X[n]  and 
X[n  —  l ]  is  described  by  the  time- dependent  covariance  function ,  which  is  defined  by 

Cx[n;  l ]  =  £{(A"[n]  —  mx[n])(X[n  —  l]  —  mx[n  —  /])*}.  (2-13) 

The  relationship  between  the  correlation  function  and  the  covariance  function  is 

Rx [n;  1}  =  C x [n;  l }  +  mx [' n]m*x [n  —  l\,  (2-14) 

hence  when  X[n]  is  a  zero-mean  random  process, 

Rx[n;  l }  =  Cx[n ;  /]. 

If  we  consider  two  random  processes,  X[n]  and  Y [n] ,  the  correlation  between 
any  two  samples  of  the  random  processes  is  described  by  the  time- dependent  cross¬ 
correlation  function ,  which  is  defined  by 

Rxrln ;  /]  =  £{X[n]Y*[n  -  l}}.  (2.15) 

An  expression  can  be  written  for  the  time- dependent  cross-covariance  function  as 

CxY[n\  /]  =  £{(X[n\  —  mx[n])(Y[n  —  l]  —  my[n  —  /])*}.  (2-16) 

The  relationship  between  the  cross-correlation  function  and  the  cross-covariance  func¬ 
tion  is 

Rxv[n\  /]  =  C.vr[n;  l\  +  mx[n]my[n  -  l],  (2.17) 
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hence  when  X[n]  and  Y[n]  are  zero-mean  random  processes, 

RxY[n',l\  =  Cxv[n;l\. 

Two  random  processes  are  called  orthogonal  if  RxY[n',l]  =  0  and  uncorrelated  if 
CxY[n;  l }  =  0. 

2.  Stationarity  of  Random  Processes 

Recall  that  a  random  process  is  wide-sense  stationary  (WSS)  if 

1.  the  mean  of  the  random  process  is  a  constant,  rrix [n]  =  rrix,  and 

2.  the  correlation  function  is  a  function  only  of  the  spacing  between  samples,  i.e., 
Rx  [n;l]  =  Rx  [l]  ■ 

and  that  two  random  processes  are  jointly  wide-sense  stationary  (JWSS)  if 

1.  they  are  each  WSS,  and 

2.  their  cross-correlation  function  is  a  function  only  of  the  spacing  between  sam¬ 
ples,  i.e.,  RXy[u\1 ]  =  Rxy[1]- 

Under  the  assumptions  of  WSS  and  JWSS,  the  mean,  correlation  and  covari¬ 
ance  functions  are  summarized  in  Table  2.4. 

3.  Matrix  Representations  of  Moments 

Using  the  vector  representation  (2.5)  for  a  random  signal,  a  number  of  impor¬ 
tant  concepts  and  properties  can  be  defined.  The  first  moment  or  mean  of  a  random 
vector  is  defined  by 

£{X[0]}  j  I"  mx[  o] 

£{X[1]}  mx[  1] 

mx  =  £{X|  =  = 

£{X[iV-l]}  mx[N-  1] 
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Mean  Function 

mx  =  £{X[n]} 

(2.18) 

(Auto) correlation  Function 

Rx[l\  =  £{X[n\X*[n~l]} 

(2.19) 

Covariance  Function 

Cx[l ]  =  £{(X[n]  -  mx)(X[n  -  l]  -  mx)*} 

(2.20) 

Interrelation 

Rx[l]  =  Cx[l]  +  \mx\2 

(2.21) 

Cross-correlation  Function 

RxY[l]  =  e.{X[n]Y*[n-l]} 

(2.22) 

Cross-covariance  Function 

CXy[1]  =  £{(X[n]  -  mx)(Y[n  -  l]  -  mY )*} 

(2.23) 

Interrelation 

Rxy[1]  =  CXy[l\  +  mxmy 

(2.24) 

Table  2.4.  Summary  of  definitions  and  relationships  for  stationary  random  processes, 
(After  [Ref.  5]). 


which  is  completely  specified  by  the  associated  mean  function  mx\p]  in  (2.10).  If  the 
random  process  is  WSS,  then  the  mean  function  is  independent  of  the  sample  index 
and  mx  is  defined  by  a  vector  of  constants 


rrix 


mx 


mx 


(2.26) 


mx 


The  correlation  matrix  represents  the  complete  set  of  second  moments  for  the 
random  vector  and  is  defined  by 


Rx  =  £{XX*T}. 


(2.27) 
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The  correlation  matrix  thus  has  the  explicit  form 


Rx  = 


£{|X[0]|2} 

£{X[1]X*[0]} 


£{X[0]X*[1]} 

£{|X[1]|2} 


Z{X[0]X*[N-  1]} 
E{X[1]X*[N-  1]} 


(2.28) 


£{|X[iV  —  1]X*[0]} 

Ry[0;0] 

Rx[  l;  l] 


£{X[N  -  1]X*[1]} 

R.y[0;—  1] 

Rx[l;0] 


.  e{\x[N-i}\2} 

Rx  [0;  —  N  +  1] 
Rx[l)~N] 


(2.29) 


RX[N  -  1;  TV  -  1]  RX[N  -  1;  N  -  2] 


Rx[N  —  1;  0] 


which  is  completely  specified  by  the  associated  correlation  function  Rx[n]  l ]  in  (2.12). 
If  the  random  process  is  WSS,  then  the  correlation  is  a  function  of  only  the  sample 
spacing  and  has  the  form  of  a  Toeplitz  matrix: 


Rx  [0] 

Rx[- 1] 

Rx[~  2]  ... 

Ry[-Y+1] 

Rx[  1] 

Rx  [0] 

Rx[- 1] 

Ry[-Y] 

Rx[  2] 

Rx[  1] 

(2.30) 

Rx[N-  1] 

Rx[N-  2] 

...  Rx[  1] 

Rx  [0] 

This  matrix  is  completely  specified  by  the  associated  correlation  function  Rx[l ]  in 
(2.19). 

The  cross-correlation  matrix  represents  the  complete  set  of  second  moments 
between  two  random  vectors  X  G  lw  and  Y  G  Mm  and  is  defined  by 


Rxy  =  £{XY't1, 


(2.31) 


and  the  associated  correlation  matrix  has  the  form 


Rxy[  0;  0]  Rxy[  0;  —1]  •  •  •  Rxy[  0;  —M  +  1] 

Rxy[  l;  l]  Rxy[  l;  0]  ...  i?.vy  [l;  —  M] 

RxY  =  , 

RYy[iV-l;lV-l]  Rxy  [IV  -  1;  N  -  2]  ...  i?Yy  [N  -  1;  IV  -  M\ 

(2.32) 

which  is  completely  specihed  by  the  associated  cross-correlation  function  /?Yy  [n;  l]  in 
(2.15).  In  general,  Rxy  is  not  a  square  matrix  (unless  N  =  M ).  If  the  associated 
random  processes  are  JWSS,  then  the  cross-correlation  is  a  function  of  only  the  sample 
spacing 

Rxy[0]  Rxy[~  1]  •••  Rxy[~M+  1] 

Rxy[  1]  Rxy  [0]  i  Rxy[~M ] 

Rxy  =  Rxy[  2]  Rxy[  1]  '•  •••  >  (2.33) 

Rxy[N-  1]  Rxy[N-  2]  ...  Rvy[iV-M]_ 
which  is  completely  specihed  by  the  associated  correlation  function  Rx[l]  in  (2.22). 
In  general,  such  matrices  will  exhibit  Toeplitz  structure  but  will  not  be  Hcrmitian 
symmetric  [Ref.  5].  Similar  expressions  and  statements  can  be  made  concerning 
the  cross-covariance  matrix  and  function.  The  essential  definitions,  properties,  and 
relations  for  the  quantities  discussed  in  this  section  are  listed  in  Table  2.5. 

4.  Reversal  of  First  and  Second  Moment  Quantities 

Since  the  operations  of  expectation  and  reversal  commute,  we  have  the  follow¬ 
ing  relations  for  the  first  and  second  moment  quantities 

mx  =  £{X}  =  mx,  (2.34) 

and 

Rx  =  £{XX*T}  =  Rx  (Cx  =  Cx).  (2.35) 
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Mean 

mx  =  £{X} 

(Auto)  correlation 

Rx  =  £{XX*T} 

Covariance 

Cx  =  £{(X-mx)(X-mx)*T} 

Interrelation 

Rx  =  Cx  +  mxmx+i 

Cross-correlation 

Rxy  =  £{XY*t} 

Cross-covariance 

Cxy  =  £{(X  —  mx)(Y  —  mY)*T} 

Interrelation 

Rxy  =  Cxy  +  mxmy*2 

Symmetry 

i'x  —  i>X  ,  ^x  —  ^x 

Relation  of  Rxy  and  Cxy 

13  _  13  *T  *T 

i'XY  —  1'  YX  ,  ^XY  ~  C/yX 

Table  2.5.  Summary  of  useful  definitions  and  relationships  for  random  processes, 
(After  [Ref.  5]). 

Further,  if  Rx  (Cx)  is  a  Toeplitz  correlation  (covariance)  matrix  corresponding  to  a 
WSS  random  process,  it  follows  that 

Rx  =  Rx  (2.36) 


D.  NUMBER  THEORY 

Number  theory,  “. . .  the  branch  of  mathematics  concerned  with  the  study  of 
the  properties  of  the  integers  [Ref.  55],”  is  a  natural  framework  for  the  analysis  of 
discrete-time  systems,  where  the  independent  variables,  by  definition,  are  integers. 
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In  particular,  since  in  this  analysis  of  multirate  systems,  notions  of  divisibility,  factor¬ 
ization  and  congruence  are  integral,  the  ensuing  discussion  is  provided  to  introduce 
and  define  these  and  related  concepts  [Ref.  55,  56,  57,  58]. 

1.  Division  Algorithm  Theorem 

The  elementary  operation  of  division  forms  the  basis  of  much  of  what  is  to 
follow  and  is  expressed  by  the  division  algorithm  theorem. 

Theorem  1.  Let  a  and  h  be  integers  with  a  >  0.  Then  there  exists  unique  integers 
q  and  r  satisfying 

b  =  qa  +  r,  0  <  r  <  a,  (2.37) 

where  q  is  called  the  quotient  and  r  is  called  the  remainder. 

The  proof  of  this  can  be  found  in  many  texts,  e.g.,  [Ref.  55,  56,  57]. 

Example  1.  A  specific  example  demonstrating  the  division  algorithm  theorem  is  pro¬ 
vided.  Given  integers  a  =  3  and  b  =  22,  we  find  unique  integers  q  =  7  and  r  —  1  that 
satisfy  (2.37).  The  quotient  is  q  =  7;  the  remainder  is  r  —  1.  ■ 

2.  Divisibility 

Definition  13.  Let  a  and  b  be  integers.  Then  a  divides  b ,  written  a\b.  if  and  only  if 
there  is  some  integer  c  such  that  b  =  ca.  When  this  condition  is  met,  the  following 
are  equivalent  statements:  (i)  a  is  a  factor  of  b ,  (ii)  b  is  divisible  by  a,  and  (iii)  b  is  a 
multiple  of  a.  If  a  does  not  divide  b,  we  write  a\b. 

Example  2.  This  example  illustrates  the  concept  of  divisibility  for  a  number  of  integer 
pairs. 

3 1 12,  7| 21,  9 1 108, 12 1 144; 

4  {5,7  {8,8  {7,3  {22.  ■ 
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a.  Greatest  Common  Divisor 

Definition  14.  Let  a  and  6  be  integers.  The  integer  d  is  called  the  greatest  common 
divisor  of  a  and  b,  denoted  by  gcd(a,  b),  if  and  only  if 

1.  d  >  0, 

2.  d\a  and  d\b, 

3.  whenever  e| a  and  e| b,  we  have  e\d. 

The  integers  a  and  b  are  said  to  be  relatively  prime  if  gcd(a,  b)  =  1. 

Example  3.  A  few  examples  demonstrating  the  greatest  common  divisor: 

If  a  =  3  and  6  =  4,  then  d  =  gcd(3, 4)  =  1  (3  and  4  are  relatively  prime), 

If  a  =  12  and  b  =  15,  then  d  =  gcd(12, 15)  =  3, 

If  a  =  25  and  b  =  55,  then  d  =  gcd(25,  55)  =  5.  ■ 

b.  Least  Common  Multiple 

Definition  15.  Let  a  and  b  be  positive  integers.  The  integer  m  is  called  the  least 
common  multiple  of  a  and  b,  denoted  by  lcm(a,  b),  if  and  only  if 

1.  m  >  0, 

2.  a\m  and  b\m,  and 

3.  if  n  is  such  that  a \ n  and  b\n,  then  m\n. 

The  least  common  multiple  can  be  expressed  as 

ah 

1cm  (a,  6)  =  — —  — .  (2.38) 

gcd(a,  b) 

Example  4.  A  few  examples  demonstrating  the  least  common  multiple: 

If  a  =  3  and  6  =  4,  then  m  =  lcm(3, 4)  =  12, 

If  a  =  12  and  6  =  15,  then  m  =  lcm(12, 15)  =  60, 

If  a  =  25  and  6  =  55,  then  m  =  lc-m(25,  55)  =  275.  ■ 
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Also  note  that  the  least  common  multiple  is  associative  and  therefore, 

lcm(a,  b,  c)  =  lcm(lcm(a,  b),c)  =  lcm(a,  lcm(6,  c).  (2.39) 

3.  Greatest  Integer  Function 

The  greatest  integer  function,  often  called  the  floor  function,  is  defined  as 
follows. 

Definition  16.  For  any  x  G  M,  the  greatest  integer  function  evaluated  at  x  returns 
the  largest  integer  less  than  or  equal  to  x.  This  is  sometimes  referred  to  as  the  integral 
part  of  x.  The  function  will  be  denoted  as  |_^J  • 

Example  5.  The  following  examples  illustrate  this  definition, 

L2.7J  =  2, 

L0.9J  =  0, 

L-0.3J  =  -1.  ■ 

Note  that  the  floor  function  satisfies  the  following  identity 

\_x  +  k\  =  [xj  +  k,  for  k  G  Z.  (2.40) 

4.  Congruence 

If  a  is  fixed  in  (2.37),  then  there  are  an  infinite  number  of  choices  of  h  for  which 
the  remainder  r  is  the  same.  In  this  context,  a  is  called  the  modulus,  the  choices  of  b 
are  said  to  be  congruent  modulo  a,  and  the  remainder  is  called  the  common  residue 
modulo  a  or  simply  the  common  residue  [Ref.  58].  This  concept  of  congruence  is 
formalized  with  the  following  definitions. 

Definition  17.  Let  n  be  a  positive  integer.  The  integers  x  and  y  are  “congruent 
modulo  n”  or  ux  is  congruent  to  y  modulo  n”,  denoted  x  =  y  (mod  n) ,  provided  that 
x  —  y  is  divisible  by  n.  If  x  and  y  are  not  congruent  modulo  n  or  x  is  not  congruent 
to  y  modulo  n,  we  write  x  ^  y  (mod  n). 
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Example  6.  We  demonstrate  the  concept  of  congruence  with  a  few  examples. 

8  =  5  (mod  3), 

14  =  2  (mod  12), 

49  =  42  (mod  7).  ■ 

Example  7.  In  the  following  example,  n  =  2,  and  there  are  two  sets  of  integers  that 
are  congruent  modulo  2,  the  even  integers  and  the  odd  integers. 

{. . . ,  —4,  —2,  0,  2, 4, . . .}  are  congruent  to  0  (mod  2), 

{. . . ,  —3,  —1, 1,  3,  5, . . .}  are  congruent  to  1  (mod  2).  ■ 

Example  8.  In  this  example,  n  —  3,  and  there  are  three  sets  of  integers  that  are 
congruent  modulo  2. 

{. . . ,  —6,  —3,  0,  3,  6, . . .}  are  congruent  to  0  (mod  3), 

{. . . ,  —5,  —2, 1, 4,  7, . . .}  are  congruent  to  1  (mod  3), 

{. . . ,  —4,  —1,  2,  5,  8, . . .}  are  congruent  to  2  (mod  3).  ■ 

Definition  18.  If  x  =  y  (mod  n),  then  y  is  called  a  residue  of  x  modulo  n.  Further¬ 
more,  if  0  <  y  <  n,  then  y  is  called  the  common  residue  of  x  modulo  n,  or  simply  the 
common  residue. 

Example  9.  Referring  to  Example  6,  we  point  out  the  associated  residues. 

5  is  a  residue  of  8  modulo  3, 

2  is  the  common  residue  of  14  modulo  12,  and 
42  is  a  residue  of  49  modulo  7.  ■ 

Definition  19.  The  set  of  integers  An  =  {0, 1, . . . ,  n  —  1}  is  called  the  set  of  “least 
positive  residues  modulo  n” 

At  times,  it  is  necessary  to  extract  the  common  residue  [Ref.  58].  This  oper¬ 
ation  is  denoted  by  (•)„  and  is  defined  as 

X 

y  =  (x)n  =  X-  -  n,  (2.41) 

n 

where  y  is  the  common  residue  of  x  modulo  n,  and  |_-J  is  the  floor  operation. 
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Example  10.  A  few  examples  of  extracting  the  common  residues  of  x  modulo  n. 


y  =  (22)3 
y  =  (14)4 


22-21 


1, 


14-  12  =  2. 


E.  CHAPTER  SUMMARY 

This  chapter  introduces  various  fundamental  signal  processing  and  mathemati¬ 
cal  concepts  required  for  theoretic  and  application-related  developments  in  subsequent 
chapters.  Further,  for  the  purposes  of  consistency,  accuracy  and  ease  of  understand¬ 
ing,  conventions  and  notation  are  also  established. 

The  taxonomy  of  signals  and  sequences,  their  various  definitions,  and  associ¬ 
ated  notations  are  presented.  Of  particular  relevance  is  the  discussion  on  discrete- 
domain  signals  and  their  sequence  representation,  which  form  the  most  basic  con¬ 
stituent  of  any  multirate  system  (Chapter  III). 

Many  concepts  from  linear  algebra  are  recalled,  including  the  concept  of  a 
random  vector  and  the  reversal  of  a  vector  or  matrix.  Further,  the  linear  algebraic 
concept  of  the  Kronecker  product  is  discussed,  which  is  useful  in  the  matrix  represen¬ 
tation  of  various  multirate  operations  in  Chapter  III  and  the  multirate  Wiener-Hopf 
equations  in  Chapter  IV.  Finally,  the  Frobenius  inner  product  is  introduced,  which 
provides  a  compact  representation  of  the  two-dimensional  linear  estimate  required  for 
image  reconstruction  (Chapter  V). 

In  the  analysis  of  random  processes,  the  second-moment  properties  are  fre¬ 
quently  used.  Since  they  are  essential  to  the  development  of  optimal  estimation 
theory,  the  analysis  and  various  definitions  and  relationships  are  reviewed  in  this 
chapter. 
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Finally,  several  topics  in  number  theory  are  presented,  which  have  great  utility 
in  developing  the  theory  of  multirate  systems  and  characterizing  the  relationships 
between  constituent  signals  and  the  related  multirate  system  (Chapter  III). 
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III.  MULTIRATE  SYSTEMS:  CONCEPTS  AND 

THEORY 

In  this  chapter,  we  develop  the  theory  of  multirate  systems,  which  establishes 
the  fundamental  relationships  in  a  multirate  system,  and  culminates  in  a  systematic 
framework  for  their  analysis.  These  results  lead  to  representation  of  the  various 
signals  in  a  multirate  system  on  a  common  domain,  system  and  impulse  response 
formulations  at  both  the  signal-  and  system-level,  linear  algebraic  representation  of 
multirate  operations,  and  ultimately,  as  presented  in  Chapter  IV,  development  of 
multirate  signal  estimation  theory. 

A.  INTRODUCTION 

In  many  digital  signal  processing  (DSP)  applications,  the  systems  involved 
must  accommodate  discrete-domain  signals  that  are  not  all  at  the  same  sampling 
rate.  For  instance,  consider  a  system  in  which  the  signals  at  the  source  and  desti¬ 
nation  have  different  sampling  rate  requirements.  An  example  of  this  occurs  when 
recording  music  from  a  compact  disc  (CD)  system  at  44.1  kHz  to  a  digital  audio  tape 
(DAT)  system  at  48  kHz.  Another  application  might  involve  systems  that  incorporate 
several  signals  collected  at  different  sampling  rates.  Sensor  networks,  many  military 
weapon  and  surveillance  systems,  and  various  controllers  process  data  from  various 
sensors,  where  the  information  from  each  sensor  might  be  collected  at  a  different  rate. 
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Further,  a  system  may  be  at  a  rate  that  is  inefficient,  and  sampling  rate  conversion 
may  be  required  to  reduce  the  rate  of  the  system  because  “oversampling”  is  wasteful 
in  terms  of  processing,  storage  and  bandwidth. 

B.  MULTIRATE  SYSTEMS 

The  various  ideas  described  in  this  chapter  follow  [Ref.  3,  59];  however,  many 
important  extensions  are  made  to  align  results  with  the  theory  of  multirate  systems  as 
developed  here.  A  multirate  system  will  be  defined  as  any  system  involving  discrete- 
domain  signals  at  different  rates.  Recall  from  Chapter  11  that  we  will  use  sequence 
notation  (i.e. ,  x[n])  and  different  index  values  (n,  m,  etc.)  to  denote  discrete-domain 
signals  at  different  rates.  Figure  3.1  depicts  a  notional  multirate  system  where  the 
input,  output  and  internal  signals  are  at  different  rates.  A  specific  example  of  a  mul- 


Figure  3.1.  Notional  multirate  system  where  input,  output,  and  internal  signals  are 
at  different  rates.  (From  [Ref.  3]). 

tirate  system  is  the  subband  coder  illustrated  in  Figure  3.2.  The  signals  x[n\  and  y[n] 
at  the  input  and  output  of  the  system  are  at  the  original  sampling  rate  while  some  of 
the  internal  signals  {y\  [mi]  and  7/2  [^2])  are  at  lower  rates  produced  through  filtering 
and  decimation. 
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x[n\ 


y[n] 


HAz) 

— 

ILj 

— ► 

Coding 

- ► 

T  Lj 

— ► 

GAz) 

y7[«h] 

Z,[n\ 

H2(z ) 

— 

il2 

— ► 

Coding 

- ► 

— ► 

G2(z) 

y2[m2. ]  Z2[n] 

Figure  3.2.  Simple  subband  coding  system. 


1.  Intrinsic  and  Derived  Rate 

The  notion  of  rate  was  introduced  in  Chapter  II  and  is  part  of  the  description 
of  any  discrete-domain  signal.  The  rate  associated  with  a  particular  signal  may  be 
a  result  of  sampling  an  analog  signal  or  a  result  of  operations  on  sequences  in  the 
system.  These  issues  are  discussed  below. 

a.  Intrinsic  Rate 

A  discrete-domain  signal  may  be  derived  from  an  analog  signal  by  pe¬ 
riodic  or  uniform  sampling  described  by 

x[n\  =  x(nTx)  =  x{t)\t=nTx  n  G  Z.  (3.1) 

Here,  x[n]  is  the  discrete-domain  sequence  obtained  by  sampling  the  analog  signal 
x(t)  every  Tx  seconds.  This  concept  is  depicted  in  Figure  3.3. 


x(t) 

— Y— 

x[n\  =  x{nTx ) 

T 

X 

Figure  3.3.  An  analog  signal  sampled  with  a  sampling  interval  of  Tx. 
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The  sampling  interval  Tx  and  its  reciprocal,  the  sampling  rate  Fx,  are  related  by 

Fx  =  yx  (3-2) 

In  this  context,  we  say  x[n]  is  at  a  rate  Fx.  The  rate  associated  with  the  sequence, 
therefore,  is  the  rate  at  which  its  underlying  analog  signal  was  sampled  and  is  referred 
to  as  its  intrinsic  rate. 

b.  Derived  Rate 

The  process  of  sampling  rate  conversion  provides  another  context  in 
considering  the  notion  of  rate  or  sampling  rate  in  multirate  systems.  The  two  basic 
operations  in  sampling  rate  conversion  are  downsampling  and  upsampling  (with  ap¬ 
propriate  filtering).  These  operations  are  depicted  by  the  blocks  shown  in  Figure  3.4, 
and  they  are  mathematically  represented  by 

y[n\  =  x[Mn\,  (3.3) 


where  n  is  an  integer,  in  the  case  of  downsampling,  and 

[x  [-1  ,  if  n\L  ; 

y[n]=l  LLJ  (3.4) 

I  0,  otherwise, 

in  the  case  of  upsampling.  Figures  3.5  and  3.6  graphically  depict  the  downsampling 
and  upsampling  operation,  respectively,  for  M  =  L  =  2. 


Figure  3.4.  Basic  operations  in  multirate  signal  processing,  downsampling  and  up- 
sampling. 
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Figure  3.5.  An  example  of  the  downsampling  operation  (3.3),  M  =  2. 


Figure  3.6.  An  example  of  the  upsampling  operation  (3.4),  L  =  2. 
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Note  that  both  operations  are  performed  exclusively  in  the  digital  domain.  The 
resulting  signals  have  no  intrinsic  rate,  but  the  rate  is  derived  from  the  rate  of  the 
input  signal.  For  downsampling,  the  output  rate  Fy  is  given  by 


F  -  — 
y  M' 


while  for  upsampling  the  rate  is  given  by 


(3.5) 


Fy  —  LFX.  (3.6) 

The  parameter  M  in  downsampling  is  called  the  decimation  factor  while  the  parame¬ 
ter  L  may  be  called  the  upsampling  factor.  Thus,  downsampling  results  in  a  reduction 
of  sampling  rate  by  a  factor  of  M,  and  upsampling  results  in  an  increase  in  sampling 
rate  by  a  factor  of  L. 

ft  will  be  seen  later  that  other  operations  more  general  than  downsam¬ 
pling  and  upsampling  can  result  in  rate  changes.  These  more  general  operations  will 
be  represented  by  linear  periodically- varying  filters  (see  Section  III.D.2).  The  out¬ 
puts  of  these  filters  have  no  intrinsic  rate  but  have  a  derived  rate  associated  with  the 
operation  that  is  performed. 


C.  CHARACTERIZATION  OF  MULTIRATE  SYSTEMS 

In  the  discussion  of  multirate  system  concepts  and  associated  theory,  it  is 
necessary  to  further  develop  terminology,  characterize  such  systems  and  develop  a 
conceptual  framework  by  which  further  analysis  and  extension  can  be  based.  In  this 
section,  the  concepts  and  terms  are  introduced. 

1.  System  Rate 

Consider  a  multirate  system  with  just  two  signals  at  different  sampling  rates, 
X\  and  X2 ■  Although  it  is  not  strictly  necessary,  the  discussion  can  be  more  easily 
motivated  if  it  is  assumed  that  each  signal  is  derived  by  sampling  an  underlying  analog 
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signal  as  shown  in  Figure  3.7.  It  will  be  assumed  that  the  sampling  rates  F\  and  F2 
are  integer-valued.  While  the  treatment  could  be  generalized  to  the  case  where  the 
rates  are  rational  numbers,  the  assumption  of  integer- values  simplifies  the  discussion 
and  is  quite  realistic  for  practical  systems.  The  corresponding  discrete-domain  signals 


h(0 

— Y— 

xT](t) 

Fi 

x,(0 

Y 

xTi  (?)  <->  x2[m2] 

F2 

Figure  3.7.  Two  signals  sampled  at  different  sampling  rates. 


and  xt2  at  the  output  of  the  samplers  are  defined  at  points  on  their  respective 
domains 


Tti  =  {nT\  ;  n  G  Z}, 

(3.7) 

and 

T t2  =  {nT2  ;n  G  Z}, 

(3.8) 

where  T\  —  —  and  T2  =  — .  The  discrete-domain  signals  are  represented  in  Figure 
F)  F-2 

3.8  as  sequences  with  different  index  values  x[n\\  and  a; [rz2]  indicating  the  different 

sampling  rates. 

Note  that  there  is  some  common  domain 

Ty  =  { nT  ;  n  G  Z}, 

(3.9) 

with  some  maximum  sampling  interval  T  in  which  the  samples  of  both  x\  and  x 2  can 
be  represented.  I11  other  words,  C  and  C  vFy. 

The  sampling  interval  T  in  (3.9)  will  be  called  the  system  sampling  interval  or  clock 
interval.  We  can  state  the  following  theorem. 
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Figure  3.8.  Two  signals  sampled  at  different  integer-valued  sampling  rates.  A  periodic 
correspondence  between  indices  can  be  observed  (as  indicated  by  the  dashed  lines). 
The  system  grid  is  represented  by  the  line  segment  at  the  bottom  of  the  figure  and 
is  derived  from  the  set  of  hidden  and  observed  samples  of  the  associated  underlying 
analog  signals.  Open  circles  represent  “hidden”  samples. 


Theorem  2.  The  system  sampling  interval  is  given  by  T 
integer,  and 

F  =  lcm(Fi,  F2). 

Here,  F  is  called  the  system  rate  or  the  fundamental  rate. 

Proof.  From  the  definition  of  the  least  common  multiple  (Definition  15): 

1.  By  definition,  F  must  be  a  positive  integer; 

2.  TTl  C  %  =>  Fi\F,  TT2  C  %  =>•  F2\F] 

3.  Since  T  is  the  maximum  sampling  interval  in  which  samples  of  both  x\  and 
x2  can  be  represented  in  ,  this  implies  that  F  =  =  is  the  related  minimum 
sampling  rate,  and  the  third  condition  of  Definition  15  is  met. 


=  =,  where  F  is  a  positive 


(3.10) 


□ 
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We  also  define  a  system  grid  as  the  representation  of  the  set  on  the  real  line, 
as  seen  in  Figure  3.8.  Note  that  the  samples  of  a  signal  corresponding  to  the  system 
grid  can  be  viewed  as  a  set  of  hidden  and  observed  samples  of  the  underlying  analog 
signal.  For  example,  for  a  given  signal  xi,  the  observed  samples  are  the  samples  of 
xi(t)  that  correspond  to  the  set  'F^  defined  in  (3.7)  while  the  hidden  samples  of  x\ 
are  the  samples  of  x\ (t)  that  correspond  to  the  set  —  'Ft/- 

It  is  frequently  useful  to  represent  all  signals  at  the  “system  level”  defined 
by  T  and  F.  The  system  level  is  referred  to  as  the  “fundamental  layer  ”  in  [Ref. 
59].  Sequences  associated  with  the  system  level  will  be  denoted  by  a  symbol  with  an 
overbar,  as  in  x\  [n]  and  x2  [ n ]  and  a  common  index  n.  This  point  is  developed  further 
in  Section  1I1.C.5. 

For  a  multirate  system  comprised  of  M  signals,  the  definition  (3.10)  can  be 
extended  (see  (2.39))  as 

F  =  lcm(Fi,  /•[>,.. . .  Fm),  (3.11) 

with  T  =  =. 

F 

2.  Decimation  Factor 

Recall,  from  Section  II. D. 2(b),  that  the  least  common  multiple  m  is  a  number, 
which  is  a  multiple  of  both  associated  integers  a  and  b,  therefore,  a\m  and  b\m. 
Further,  from  Definition  13,  the  condition  a\m  implies  that  m  =  cia,  and  likewise 
b\m  implies  that  m  =  C2 b,  where  c\  and  C2  are  constants. 

We  can  apply  these  results  to  multirate  systems.  Since  F  is  the  least  common 
multiple  of  F\  and  F2,  it  follows  that  there  exists  integer  constants  K\  and  K2  such 
that 

F  =  K  i  F\  and  F  =  K2  F2 .  (3-12) 

1The  “difference”  A  —  B  of  two  sets  A  and  B ,  where  B  C  A,  is  defined  as  A  —  B  =  {x  €  A\x  (j  B}. 
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Rearranging  Equation  (3.12)  yields  an  expression  for  these  constants  as 

Kx  =  and  K2  =  ^.  (3.13) 

e  i  r2 

Further,  for  a  multirate  system  comprised  of  M  signals,  (3.13)  can  be  extended  as 


F 

Ki  =  —  where  i  —  1, . . . ,  M. 
Fi 

Notice  that  if  (3.14)  is  rearranged  to  the  form  of  (3.5): 

F 


(3.14) 


F,;  = 


AV 


we  explicity  see  that  the  system  sampling  rate  F  is  reduced  by  a  factor  of  A"*.  There¬ 
fore,  Ki  is  defined  as  the  system  decimation  factor,  or  simply  the  decimation  factor, 
for  the  ith  signal. 

As  a  consequence,  the  following  observation  can  be  made  and  extended  to 
any  number  of  signals.  If  the  underlying  analog  signals  X\  and  x2  are  sampled  at 
the  system  rate  F  and  then  are  decimated  by  their  respective  decimation  factors  K± 
and  K2 ,  the  resultant  discrete-domain  signals  are  the  signals  obtained  by  sampling 
xi  and  x2  at  the  original  rates  F\  and  F2,  respectively  (see  Figure  3.9).  In  terms  of 
the  hidden  and  observed  samples  of  x\  and  x2,  the  decimation  factors  offer  a  way  to 
relate  the  set  of  observed  samples  to  the  set  of  all  samples  (hidden  and  observed)  at 
the  system  level. 

If  (3.14)  is  expressed  in  terms  of  sampling  intervals,  we  have 

(3.15) 

where  the  decimation  factor  is  the  ratio  of  the  duration  of  the  ith  signal’s  sample 
interval  to  the  duration  of  its  associated  clock  interval.  Thus,  the  decimation  factor 
Ki  can  also  be  viewed  as  the  number  of  clock  intervals  or  system  samples  between 
samples  of  the  ith  signal  when  sampled  at  A). 


T 

Ki  =  =, 

T 


3.  System  Period 

Again,  consider  a  multirate  system  comprised  of  two  signals  x\  and  x2  sampled 
at  integer- valued  sampling  rates  F\  and  F2,  respectively.  If  samples  of  these  signals 
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Figure  3.9.  Signals  sampled  at  the  system  rate  and  decimated  by  their  respective 
decimation  factors  yield  the  original  discrete-domain  signals. 

align  on  the  related  system  grid  at  some  time  t  =  to,  then  their  samples  realign 
at  some  later  time  t  —  t\  —  to  +  T  and  at  regular  intervals  thereafter  (see  Figure 
3.10).  The  smallest  positive  value  of  T  for  which  the  realignment  occurs  is  called 
the  system  period  (not  to  be  confused  with  the  system  sampling  interval  T). Notice 
that  T  =  t\  —  to  and  if  t\  and  to  are  expressed  in  terms  of  their  discrete-domain 
representations,  we  can  write  T  =  (jii  —  no)T;  thus,  we  define 

N  =  n i  -  no, 

which  is  called  the  discrete  system  period.  Observe  that  N  represents  the  number 
of  system  samples  between  sample  realignments,  ft  relates  to  the  system  period  (in 
seconds)  as 

T  =  NT.  (3.16) 


Now,  we  define  Mi  to  be  the  number  of  signal  samples  per  period.  Therefore, 
we  relate  Mi  and  M2  to  the  discrete  system  period  by  M\T\  =  NT  =  M2X2.  Recalling 
(3.14),  we  write 

N  =  M  i  K]  =  M2K2 ,  (3.17) 

and  recalling  Definition  (II.G.b.15),  we  can  write 

N  =  lcm(/l1,  K2),  (3.18) 
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*i  (0 


Figure  3.10.  Two  signals  sampled  at  different  integer- valued  sampling  rates.  Observe 
the  periodic  alignment  between  indices,  (After  [Ref.  3]). 


an  explicit  expression  for  the  discrete  system  period.  Also  notice  from  (3.17)  that 


Mi 


—  and  M2 
Ai 


N 

AV 


(3.19) 


For  a  multirate  system  comprised  of  M  signals,  the  definition  of  (3.18)  can  be 
extended  (see  (2.39))  as 


N  =  lcm(Ad,  K2, ,  Km)- 


(3.20) 


4.  Maximally-decimated  Signal  Set 

Consider  a  discrete-domain  signal  y[n)  at  rate  Fy.  If  y[n]  is  downsampled  by 

a  factor  of  M,  after  M  —  1  successive  translations,  a  set  of  related  discrete-domain 

signals  results,  designated  {xo[m],  £i[m], . . . ,  sm-iH).  Note  that  the  constituent 

p 

discrete-domain  signals  are  at  rate  Fx  =  An  example  is  shown  in  Figure  3.11  for 
M  =  3.  If  the  discrete-domain  signals  associated  with  y[n]  are  described  by 

Xi[m\  =  y[i  +  nM]  i  —  0, 1, . . . ,  M  —  1,  (3.21) 
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Figure  3.11.  Example  of  a  3-fold  maximally  decimated  signal  set. 


then  the  resultant  signals  will  be  called  a  maximally- decimated  signal  set  with  down- 
sampling  factor  M  or  an  M-fold  maximally -decimated  signal  set. 

5.  Representation  of  Signals  in  Multirate  Systems 

In  the  analysis  of  multirate  systems,  there  is  a  need  to  relate  signals  to  a 
common  scale.  This  scale  is  represented  by  the  system  grid  and  is  discussed  in  terms 
of  the  Tyr  domain,  where  T  is  the  system  sampling  interval.  The  sample  indices  of  a 
signal  sampled  at  the  system  rate  correspond  to  the  integer  multiples  of  T .  In  this 
context,  every  signal  in  a  multirate  system  can  be  represented  on  the  system  grid. 

Consider  a  multirate  system,  with  a  system  sampling  interval  T,  containing  a 
particular  discrete-domain  signal  Xtx  ,  among  others.  This  signal  can  be  represented 
by  the  sequence 

x[mx]  =  xTx(mxTx), 

where  xtx  is  the  discrete-domain  signal  and  Tx  is  the  associated  sampling  interval. 
Note  that  the  sequence  and  its  sampling  rate  is  defined  by  the  discrete-domain  signal 
xTx{t)  for  t  =  mxTx  and  not  by  the  analog  signal  x(t),  which  may  not  be  directly 


49 


observable  or  may  not  even  exist!  We  refer  to  this  as  a  signal-rate  or  signal  level 
representation  of  x  and,  in  this  case,  say  that  x  is  at  its  native  rate.  Now  define  x y 
to  be  the  associated  system-rate  or  system-level  representation  of  x.  This  signal  can 
be  represented  by  the  sequence 

x[n]  =  XTp{nT). 

We  can  relate  the  two  representations  by  recalling  from  (3.15)  that  Tx  =  TKX  and 
noting  that 

xTx  (mxTx)  =  xT 'x{mxTKx)  =  xT(mxKxT )  =x[Kxmx\. 

Therefore,  a  discrete-domain  signal  at  its  native  rate  can  be  represented  at  the  system 
rate  or  system  level  by 

x[mx]  =  x[Kxmx],  (3.22) 

where  Kx  is  the  decimation  factor  of  signal  x.  These  concepts  are  illustrated  in 
Figure  3.12  for  a  multirate  system  comprised  of  two  signals  x  and  y.  The  native  rate 
corresponds  to  a  signal’s  “original”  sampling  rate,  and  the  system  rate  corresponds 
to  the  system  rate  determined  from  (3.10). 

6.  Summary  of  Multirate  Relationships 

For  convenience,  the  representations  and  relationships  developed  in  this  chap¬ 
ter  regarding  multirate  systems  are  summarized  in  the  following  tables.  In  Table 
3.1,  the  various  signal  representations,  their  notations  and  their  associated  rates  are 
indicated. 


x(t)  Analog 

x[n\  System-level,  system  rate,  F 

x[mx]  Signal-level,  native  rate,  Fx 

Table  3.1.  Signal  representations  in  multirate  systems. 
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Analog 


Discrete-Time 


x(t) 


y(t) 


x[mj  =  xT  ( mxTx )  Native-rate 


x[n]  =  x(Tn ) 

System-rate 

y[n\  =  y(fn) 


y[my  ]  =  yT  (myTy  )  Native-rate 


Figure  3.12.  Signal  representations  and  sampling  levels  in  a  multirate  system. 


In  Table  3.2,  the  characterization  of  and  relationship  between  signals  in  multirate 
systems  is  summarized. 


Name 

Relationship 

System  Rate 

F  =  lcm(Fi,  F2,.. 

■ ,  Fm) 

System  Sampling  Interval 

1 

T  _ 

Decimation  Factor 

Ki  =  £  =  2 

Ft  T 

System  Sample  Period 

N  —  lcm(iF1,  K2, . 

•  • ,Km ) 

Samples/Period 

System  Period 

T  =  NT  =  MiF 

Table  3.2.  Summary  of  various  relationships  pertaining  to  a  multirate  system  (M 
signals) . 
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Finally,  in  Table  3.3,  these  characterizations  and  relationships  are  summarized  in 
relation  to  signal  and  system  level  representations. 


Sampling 

Sampling 

Decimation 

Samples 

System 

Rate 

Interval 

Factor 

per  Period 

Period 

Signal 

Fx 

Tx 

Kx 

Mx 

T  =  MXTX 

System 

F 

T 

1 

N 

T  =  NT 

Table  3.3.  Parameters  pertaining  to  a  multirate  system,  (After  [Ref.  3]). 


D.  MULTIRATE  SYSTEM  THEORY 

1.  Description  of  Systems 

A  signal  processing  system  represents  the  process  for  transforming  a  signal 
into  another  signal.  The  concept  is  illustrated  in  Figure  3.13(a),  where  a  signal  x{t ) 
is  transformed  into  a  signal  y{t).  The  input  and  output  may  be  of  any  of  the  signal  or 
sequence  types  discussed  in  Chapter  II  and  need  not  be  of  the  same  type.  For  instance, 
the  input  may  be  a  discrete-domain  signal  and  the  output  may  be  an  analog  signal. 
Our  primary  concern,  however,  is  the  case  where  both  the  input  and  the  output  are 
discrete-domain  signals,  not  necessarily  with  the  same  sampling  intervals. 

Since  discrete-domain  signals  are  represented  by  sequences  together  with  a 
(known  or  implied)  sampling  rate,  it  is  appropriate  to  consider  the  properties  of 
“systems”  that  transform  sequences.  These  will  be  referred  to  as  discrete  systems 
as  shown  in  Figure  3.13(b),  and  the  transformation  will  be  represented  by  In 
addition  to  performing  a  specified  operation  on  the  input  sequence,  a  discrete  system 
provides  a  means  of  determining  the  output  sampling  rate  from  the  input  sampling 
rate.  When  input  and  output  sampling  rates  are  not  the  same,  it  is  our  convention 
to  use  different  letters  for  the  index  variable  of  the  sequence. 
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y(t) 

-► 


x(t) 


> 


T{) 


(a) 


x[n] 


y[n\ 

-► 


(b) 

Figure  3.13.  (a)  Block-diagram  representation  of  a  signal  processing  system;  (b) 

Block-diagram  representation  of  a  discrete  system. 

In  general,  a  discrete-domain  system  is  represented  by 

y [m]=7D{x.[n]}  or  7D{x.[n]}  ^  y[m],  (3.23) 

where  7d  is  a  suitable  mathematical  operator  and  x[rt]  and  y[m]  represent  the  input 
and  output  sequences  for  such  a  system.  Note  that  the  sequences  in  general  may  be 
at  different  rates. 

2.  Classification  of  Discrete  Systems 

Discrete  systems  can  be  classified  by  certain  restrictions  or  characteristics 
placed  upon,  or  observed  concerning  them.  A  number  of  such  classifications  are  impor¬ 
tant  in  signal  processing,  including:  linear/non-linear,  time-invariant /time- variant, 
causal/non-causal,  stable/unstable,  invertible/non-invertible  systems,  and  systems 
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with  memory/without  memory.  In  this  section,  we  deal  with  the  first  three  of  these 
classifications.  Others  are  not  so  important  to  our  discussion  and  conform  to  defini¬ 
tions  used  in  many  text  books  [Ref.  60,  61]. 

a.  Linearity 

A  discrete  linear  system,  is  any  system  that  satisfies  the  generalized 
superposition  principle.  If  the  input  to  such  a  system  consists  of  a  weighted  sum  of 
sequences,  then  its  output  consists  of  the  weighted  sum  of  the  system  responses  to 
each  individual  input  sequence  (superposition  of  responses).  It  is  sufficient  to  define 
linearity  in  terms  of  the  system’s  response  to  a  weighted  sum  of  two  sequences. 

Definition  20.  Let  x \  [n]  and  x2[n\  be  two  sequences  at  the  same  sampling  rate.  A 
discrete  system  is  linear  if  and  only  if 

7D{oiiXi[n]  +  a2x2  [n]}  =  a17D{x1[n]}  +  a27 D{x2[n]} ,  (3.24) 

for  arbitrary  input  sequences  xi[n]andx2[n]  and  constants  aianda2. 

Note  that,  while  the  output  of  a  linear  system  may  be  at  a  different 
rate  than  the  input,  we  make  no  attempt  to  define  linearity  in  terms  of  two  sequences 
at  different  rates.  Indeed,  the  sum  of  two  sequences  at  different  rates  has  no  obvious 
or  unique  intuitive  meaning  and  is  not  necessary  for  our  theory  of  multirate  systems. 

Any  discrete  system  that  does  not  satisfy  Definition  20  is  a  nonlinear 
discrete  system  (e.g.,  square  law  system  defined  by  y[n]  =  (x[n])2).  Time  dependence 
is  not  an  issue.  Both  downsamplers  and  upsamplcrs  (see  Section  III. B.  1(b))  are  linear 
systems. 

b.  Shift-invariance 

Systems  can  also  be  characterized  by  the  variation  in  their  input-output 
characteristics  as  time  evolves,  and  can  be  subdivided  into  shift-invariant  and  shift- 
dependent  systems  (frequently  called  time-invariant  and  time-dependent). 
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Definition  21.  Let  7d  be  a  discrete  system  such  that  the  input  and  output  sequences 
are  at  the  same  rate,  i.e.,  7o{x[n]}  =>■  {y[n]}.  Then  the  system  is  shift-invariant  if 
and  only  if 

7o{x[n  —  iV] }  y[n  —  IV],  for  all  integers  N.  (3.25) 

If  the  system  satisfies  (3.25)  only  for  a  particular  value  N  and  integer 
multiples  thereof,  the  system  is  said  to  be  periodic  with  period  N.  For  systems  with 
different  input  and  output  rates,  shift-invariance  is  not  defined.  Periodicity,  however, 
can  be  generalized  to  include  these  systems,  as  shown  below. 

c.  Periodic  Shift-invariance 

For  discrete  systems  with  input  and  output  signals  at  different  rates, 
i.e.,  T£){x[n]}  =>■  y [m],  we  define  a  particular  type  of  shift-invariance  called  peri¬ 
odic  shift-invariance  (PSI).  A  system  that  observes  this  property  is  called  a  discrete 
periodically  shift-invariant  system.  The  definition  of  this  property  follows. 

Definition  22.  Let  7d  be  a  discrete  system  such  that  if  7o{x[mf\}  =>■  {y[my\},  then 
the  system  is  periodically  shift-invariant  if  there  exists  integers  Mx  and  My  such  that 

7D{x[mx  -  Mx]}  =>  y[my  -  My\.  (3.26) 

Note  that  when  the  input  and  output  are  at  the  same  rate,  a  periodic 
system  satisfies  this  definition  with  Mx  =  My  and  a  shift-invariant  system  satisfies 
this  definition  for  all  Mx  such  that  Mx  =  My.  The  need  for  the  more  general  definition 
given  in  (3.26)  will  become  clear  further  in  our  development. 

d.  Causality 

In  defining  causality,  it  is  important  to  know  the  rates  or  sampling 
intervals  of  the  sequences  involved.  In  fact,  as  pointed  out  in  [Ref.  62],  “Causality  is 
intrinsically  related  to  the  ordering  in  time  [or  domain]  of  input  and  output  signals 
of  [a]  system.”  Ordering  is  clear  in  a  single-rate  system  as  identical  sample  indices 
in  each  associated  sequence  correspond  to  identical  points  in  a  common  domain.  In 
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multirate  systems,  sequence  sample  indices  must  be  referred  back  to  their  absolute 
scales  (i.e. ,  the  discrete-domain  signal)  to  discuss  ordering. 

Definition  23.  A  discrete  system  7d  is  causal  if  and  only  if 


y[my\  depends  only  on  x[mx  —  k]  for  k  —  {0, 1, . . .  }, 


(3.27) 


where  mr  — 


Ky'rriy 

Kr 


,  and  Kx,Ky  are  the  associated  decimation  factors. 


A  system  is  noncausal  if  it  does  not  satisfy  this  definition. 

It  can  be  seen  that  a  discrete  system  is  causal  if  and  only  if  the  discrete- 
domain  signal  yx y  ( t )  at  t  =  myTy  depends  only  on  values  of  the  discrete-domain  input 
signal  xj x{t)  for  t  =  mxTx  <  myTy.  This  set  of  input  values  is  known  as  the  region 
of  support  of  the  system.  The  concept  of  causality  is  illustrated  in  Figure  3.14.  The 
discrete  multirate  system  depicted  has  a  discrete-domain  input  signal  represented  by 
the  sequence  x[rnx]  and  an  output  signal  represented  by  the  sequence  y[my\.  For  a 


given  index  in  the  output  sequence  rnyo ,  Definition  23  requires  that  mXo  < 
for  causality. 


Ky‘my0 

Kr 


x[mf\ 


%{■} 


y[my\ 


Figure  3.14.  Concept  of  causality  in  a  discrete  multirate  system  comprised  of  a 
discrete-domain  input  signal  x[mx\  and  output  signal  y [rny] . 
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3.  Representation  of  Discrete  Linear  Systems 

Given  a  discrete  linear  system  with  input  x[m\  and  output  y[n],  with,  possibly, 
different  sampling  rates,  the  input-output  relation  or  system  response  can  be  written 
as 

OO 

y[n\  =  E  g[n,m]x[m].  (3.28) 

m=— oo 

The  term  g[n,m],  called  the  kernel  or  Green’s  function,  is  the  response  of  the  system 
at  point  n  in  the  output  sequence  to  a  unit  impulse2  applied  at  point  m  in  the  input 
sequence.  This  formulation  is  a  discrete-domain  version  of  the  continuous  model  in 
[Ref.  63,  64,  62],  and  the  corresponding  kernel  has  been  referred  to  as  the  Green’s 
function  weighting  pattern  response  [Ref.  65]. 

a.  Single-rate  Systems 

When  the  input  and  output  rates  are  the  same,  (3.28)  can  be  written 
as 

OO 

y[n\  =  E  h [n;  m\ x[n  —  m\,  (3.29) 

m=— oo 

where 

h[n;  m\  =  g[n,  n  —  m],  (3.30) 


is  called  the  shift- dependent  impulse  response.  The  system  is  causal  (see  Section  III. D. 2(d)) 
if 


h[n\  m ]  =  0  for  m  <  0. 


(3.31) 


The  system  is  periodic  (see  Definition  21)  if  there  exists  N  such  that 


h[n\  m\  =  h[n  +  N\m\  for  all  n. 


2The  unit  impulse,  also  know  as  the  Kronecker  delta  function,  is  defined  as 


<5[n] 


1,  n  =  0; 
0,  n  0. 


(3.32) 
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If  in  addition,  the  system  is  shift-invariant  (see  Definition  21),  then  the 
impulse  response  is  necessarily  independent  of  the  first  argument  (n)  and  the  output 
can  be  written  as  the  familiar  convolution  summation 

OO  OO 

y[n]  —  h[m\x[n  —  m}=  h[n  —  m]x[m].  (3.33) 

m=— oo  m=— oo 


The  system  represented  in  (3.33)  is  commonly  referred  to  as  a  filter  because  of  its 
interpretation  in  the  Fourier  domain,  although  the  term  “filter”  is  often  used  to 
apply  to  any  system;  linear  or  nonlinear,  shift-variant  or  shift-invariant,  in  the  signal 
processing  literature.  A  filter  is  said  to  be  a  finite  impulse  response  (FIR)  filter  (or 
system)  if  the  sequence  h[n]  has  finite  support  and  an  infinite  impulse  response  (HR) 
filter  (or  system),  otherwise.  The  region  of  support  or  the  support  of  a  given  sequence 
is  the  set  of  values  over  which  the  sequence  is  non-zero  [Ref.  66] . 
b.  Multirate  Systems 

In  developing  results  for  multirate  systems,  it  is  convenient  to  first 
describe  the  multirate  system  at  the  “system  level”  discussed  in  Section  III.C.5  of 
this  chapter.  When  represented  at  the  system  level,  we  can  apply  known  results  from 
the  analysis  of  single-rate  systems  to  the  multirate  system  and  then  express  those 
results  at  the  signal  level  pertaining  to  the  signals  of  interest. 

(1)  System-level  Representation  Recall  the  system  response 


equation  (3.28) 


y[mv\  =  y 'g[my,mx\x[mx\, 

mx 


(3.34) 


where  y  and  x  are  at  different  rates  represented  by  my  and  mx,  respectively.  Also, 
recall  (3.22)  that  a  discrete-domain  signal  XTx(t)  can  be  represented  at  the  system 
level 


x[mx]  =  x[Kxm  . 


Now  define  g,  the  system  kernel  or  system  Green’s  function ,  such  that 


g[my,mx]  =  g[Kymy,  Kxmx], 
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where  Kx  and  Ky  are  the  decimation  factors  for  x  and  y,  respectively.  At  other  values 
of  it  arguments,  g  is  defined  to  be  0.  Now,  (3.34)  can  be  written  as 

y[Kymy\  =  ^2  g[Kymy ,  Kxmx]  x[Kxmx\.  (3.35) 

mx 

If  g  satisfies  Definition  22  for  periodic  shift-invariance,  then  a  necessary 
and  sufficient  condition  is  that 


g \Kyrny ,  Aym^,]  g^Ky(^my  T  A/y),  Kx(^mx  -t-  Af^,)]. 

where  Mx  and  My  are  the  number  of  samples  per  system  period  for  x  and  y ,  respec¬ 
tively.  Then 


g[Ky(my  +  My ),  Kx(mx  +  Mx)\  =  g[Kymy  +  N,  Kxmx  +  N]  (3.36) 

where  we  recall  (3.17)  that  the  system  sample  period  is  given  by  N  =  KyMy  =  KXMX ; 
therefore,  ~g  has  period  N  in  both  arguments. 

Now,  we  define 

^[Ay'/rZy,  A 3,77X3,]  =  g[Kymy,  Kymy  Aym^,],  (3.37) 

or,  equivalently, 

g[Kymy,  Kxmx\  =  h[Kymy ,  Kymy  -  Kxmx], 

where  h  is  called  the  system  impulse  response.  Notice  that  since  g  satisfies  (3.36), 
then 

h[n;  l ]  =  g[n ,  n  —  l\  =  g[n  +  N,n  —  l  +  N]  =  h[n  +  IV;  l] , 

thus  h[n;  l }  is  periodic  in  its  first  argument  only.  Note,  also,  that  if  the  system  is 
causal,  then  h[n;  l]  must  be  0  for  l  <  0.  Now  we  can  write 

g[Kymy ;  Kxmx]  =  h[Kymy  +  N]  Kymy  -  Kxmx\ ,  (3.38) 

where  N  is  the  system  sample  period. 
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Now,  (3.35)  can  be  written  as 

y[Kymy\  =  22  h[Ky'my]  Kyrny  -  Kxmx\  x[Kxmx\.  (3.39) 

mx 

Since  h  is  periodic  in  N,  we  can  define 

h^\l]  =  h[n ;  l]  where  p  =  {n)N,  (3.40) 


and  N  is  the  system  sample  period.  Substituting  into  (3.39)  yields  the  general  filter 
equation 

1 22  ^(p)  ln  -  Kxmx\  x[Kxmx } ,  Ky  I n, 


y[n]  = 


i 


(3.41) 


undefined,  otherwise 

v 

where  p  =  (n)N. 

If  we  desire  a  causal  filter  (see  Section  III. D. 2(d)),  then 

KyUly 

mxTx  <  rriyTy  or  mx  <  — , 


(3.42) 


and 


(  I  KVmV  I 

L  Kx  J 

22  tn  _  Kx'mx\  x[Kxmx ] , 

mx=— oo 


Ky\n, 


(3.43) 


undefined,  otherwise 

where  p  =  ( n)N .  The  upper  limit  on  the  sum  is  insured  if  h  [l]  =  0  for  l  <  0.  Notice 
that  when  Tx  =  Ty  —>  Kx  =  Ky  (single-rate  system),  (3.41)  and  (3.43)  simplify  to  the 
expected  convolution  sums. 


Example  11.  Consider  Figure  3.15,  where  two  signals  are  represented  on  their  as¬ 
sociated  grids.  The  upper  grid  (a)  represents  a  signal  y  that  has  an  associated  dec¬ 
imation  factor  Ky  =  3,  and  the  lower  grid  (b)  represents  a  signal  x  that  has  an 
associated  decimation  factor  Kx  =  2.  Recall  (3.18)  that  the  system  period  is  given  by 
N  =  lcm (Kx,Ky);  therefore,  N  =  6. 
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y[n] 
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M  =  2  ' 
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np, 


11  12 
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(b) 


nT 
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(c) 

Figure  3.15.  (a)  Discrete-time  signal  y[n]  with  decimation  factor  Ky  =  3;  (h) 

Discrete-time  signal  x[n)  with  decimation  factor  Kx  =2;  (c)  System  grid. 


If  we  desire  the  estimate  for  y  at  my  =  4,  using  a  causal  filter,  where 
the  filter  order  P  =  3,  then  we  can  write  (3.43)  as 


LfJ=6_ 

y [12]  =  //°')[12  —  2mx ]  x[2 mx],  where  p  =  (12) 6  =  0. 

mx= 3 

therefore 

y[  12]  =  7i(0) [4]x[8]  +7^(0)[2]^[10]  +^(0)[0]  x[12]. 


Notice  that  the  linear  combination  is  in  terms  of  the  system-rate  parameters  (n  is  the 
system  time  index).  Once  y{  12]  is  computed,  recall  that  y[my\  =  y[myKy\;  therefore, 
y[4]=y[12]. 
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If  we  desire  the  estimate  for  y  at  rny  =  5, 


LfJ=7_ 

y[  15]  =  ^  h^\  15  —  2mx]x[2mx\,  where  p  =  (15) 6  =  3. 

mx=  5 

therefore 

y[  15]  =  h(3)  [5]x[10]  +  h(3)  [3]x[12]  +  h(3)  [l]x[14] .  ■ 

(2)  Signal-level  Representation  Often,  it  is  desirable  or  more 
convenient  to  deal  in  terms  of  the  actual  signal  parameters  rather  than  system  pa¬ 
rameters;  therefore,  we  develop  (3.41)  and  (3.43)  in  terms  of  the  individual  signal 
parameters.  To  do  so,  we  introduce  the  following  substitutions  and  change  of  vari¬ 
ables. 

First  let  us  examine  the  system  impulse  response  h, 
h  Kpm^  h  [Kyniy  T  TV ,  Kymy  h.r o i  r\ . 

Consider  the  first  argument,  where  Ky  has  been  factored  out  and  we  recall  (3.17) 
Kymy  +  N  =  Ky(my  +  N/Ky)  =  Ky(my  +  My). 

Since  this  argument  is  periodic, 


Ky{rny  +  N/Ky )  >  K y  (m  y)  j 


Now,  consider  the  second  argument,  where  Kx  has  been  factored  out 

"  K„rri, 


IS-yTYly  Kxmx  Kx 


'■y"°y 


Kr 


—  m. 


Recalling  that 


and 


y[my\  =  y[myKy\,  x[mx]  =  x[mxKx], 
h[Kymy]  Kxmx\  =  h[my-,mx\, 
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we  can  write  the  impulse  response  as 


h® 


KyUly 

Kx 


and  we  can  write 


where  l  =  (my)M  , 


y[™y\  =  5^  h(l) 


KyUly 

Kr 


—  m7 


x[mx\,  where  l  =  ( my ) 


MW 


If  we  introduce  the  following  change  of  variables 


k  = 


KyUly 

Kr 


Ul7 


(3.44) 


we  can  write  the  general  signal-rate  filtering  equation  as 


y\my]  =  ^l)[mx]x 

mx 

KyUly 

[  Kx  \ 

-  mx 

,  where  1=  {my)My. 

and  the  signal-rate  causal  filtering  equation  as 


OO 

y\my]  =  ^2  hWlmx\x 

KyUly 

Kr 

-  mx 

,  where  1=  {mv)My. 

mx= 0 

^  X  _ 

- 

(3.45) 


(3.46) 


The  lower  limit  (mx  =  0)  on  the  summation  is  equivalent  to  the  causal  condition 
h^[mx\  =  0  for  mx  <  0.  Notice  that  My ,  the  number  of  signal  samples  ( y )  per 
system  period,  is  the  period  of  the  filters  required  to  form  y[uiy],  referred  to  as  the 
cy clo stationary  period  in  [Ref.  9].  Additionally,  note  that  when  Tx  =  Tyi  then 
Kx  =  Ky  and  (3.45)  and  (3.46)  simplify  to  the  expected  convolution  sums  for  a 
single-rate  system. 


Example  12.  Let  us  illustrate  with  the  same  example  previously  discussed  and  illus¬ 
trated  in  Figure  3.15.  Recall  that  signals  y  and  x  have  associated  decimation  factors 

Ky  =  3  and  Kx  =  2,  respectively,  and  that  the  system  period  N  =  6.  Further,  recall 

,  „  ,  N  6 

that  My  =  —  —  =  2. 
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If  we  desire  the  estimate  for  y  at  rny  =  4,  using  a  causal  FIR  filter, 
where  the  filter  order  P  =  3,  then  we  can  write  (3.45)  as 

2 

y[ 4]  =  ^  h^[mx\x[ 6  —  mx\,  where  l  =  (4)2  =  0, 

mx=  0 

therefore 

y[  4]  =  fr(°>[0]x[6]  +  /i(0)[l]a;[5]  +  h^[2]x[4]. 

If  we  desire  the  estimate  for  y  at  rny  =  5, 

2 

y[ 5]  =  /?.<-1')[mx]x[7  —  mx],  where  l  =  (5)2  =  1, 

mx=0 

therefore 

y[  5]  =  /i^[0]x[7]  +  /i^[l]x[6]  +  [2]x[5].  ■ 

E.  MATRIX  REPRESENTATION 

In  the  analysis  of  multirate  systems,  linear  algebra  concepts  provide  a  useful 
framework  to  represent  or  model  basic  multirate  operations,  such  as  downsampling, 
upsampling  and  linear  filtering.  This  section  develops  a  systematic  methodology, 
which  is  used  in  other  sections  to  further  represent  multirate  systems,  an  extension 
of  work  presented  by  [Ref.  67]. 

1.  Decimation 

Decimation  is  the  process  of  digitally  reducing  the  sampling  rate  of  a  signal 
[Ref.  65],  and  the  basic  element  used  in  digital  systems  is  the  downsampler  or  dec- 
irnator,  shown  in  Figure  3.16.  This  element  extracts  every  Mth  sample  of  its  input 
and,  as  a  result,  decreases  the  sampling  rate  by  jj.  Its  operation  can  be  expressed 
mathematically  as 

y[n]  =  x[Mn]  (3.47) 
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x[n\ 


Figure  3.16.  M-fold  downsampler. 


where  n  is  an  integer. 

If  the  input  and  output  signals  are  represented  as  vectors,  then  this  operation 
can  be  expressed  as  a  linear  transformation  of  the  form 

y  =  7{M,  N,  x}  (3.48) 

where  T{-}  represents  an  arbitrary  linear  operator  that  is  dependent  on  the  decimation 
factor  M  and  the  order  of  its  associated  vectors.  If  the  downsampler  takes  input  vector 
x  G  RNM  where 

x  =  [  x[0]  x[l]  ...  x[M]  x[M  +  1]  ...  x[NM  —  1]  ]T 

and  produces  an  output  vector  y  G  M.N  where 

y  =  [  x[0]  x[M\  x[2M]  . . .  x[(N  -  1  )M\  }T  , 

then  the  operator  can  be  expressed  as 

7  =  Dm,  (3.49) 

and  the  transformation  can  be  written  as 

y  =  Dmx.  (3.50) 

The  constant  matrix  Dm  G  M NxNM  js  called  a  decimation  matrix  and  is  defined  in 
terms  of  a  Kronecker  product  (introduced  in  Section  II. B. 2)  as 

Dm  =  I  <8>  tT  (3.51) 

where  I  is  the  N  x  N  identity  matrix  and  ik  is  an  M  x  1  index  vector  with  a  1  in  the 
first  position  and  0’s  elsewhere. 
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As  an  example,  consider  the  case  where  M  —  3  and  x  e  M12.  The  decimation 
matrix  D3  G  M4x12  is  given  by 


D3 


1  0  0 
0  0  0 
0  0  0 
0  0  0 


0  0  0 
1  0  0 
0  0  0 
0  0  0 


0  0  0 
0  0  0 
1  0  0 
0  0  0 


0  0  0 
0  0  0 
0  0  0 
1  0  0 


which  results  in  y  6f4. 


2.  Expansion 

Expansion  and  decimation  are  dual  operations.  Expansion  is  the  process  of 
digitally  increasing  the  sampling  rate  of  a  signal,  and  the  basic  element  used  is  the 
upsampler  or  expander,  shown  in  Fig  3.17.  This  element  inserts  L  —  1  zeros  after  every 
sample  of  the  input  and,  as  a  result,  increases  the  sampling  rate  by  L.  Mathematically 
the  process  is  described  by 


y[n\ 


if  n  is  an  integer  multiple  of  L; 
otherwise, 


(3.52) 


or 

OO 

y[n]  =  E  x{k)5{n—kL).  (3.53) 

k= — 00 


x[n\ 


Figure  3.17.  L-fold  expander. 

If  the  input  and  output  signals  are  represented  as  vectors,  then  this  operation 
can  be  expressed  as  a  linear  transformation  of  the  form 


y  =  T{L,iV,x} 


(3.54) 
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where  T{-}  represents  an  arbitrary  linear  operator  that  is  dependent  on  the  upsam¬ 
pling  factor  L  and  the  order  of  its  associated  vectors.  If  the  expander  takes  input 
vector  x  e  lw  where 


=  x[0]  x[l]  . . .  x[N  —  1]  ]J 


and  produces  an  output  vector  y  e  M.^NL^  where 

L—l  L—l 

y  =  |  x[0]  0  ^  0  x[l]  0  ^  0  x[2] 


L- 1 

x[(N  -  1)]  O'...  0 


then  the  operator  can  be  expressed  as 


T=Ul,  (3.55) 

and  the  transformation  can  be  written  as 

y  =  ULx.  (3.56) 

The  constant  matrix  Ul  €  ]^NLxN  is  called  an  expansion  matrix  and  is  defined  in 
terms  of  a  Kronecker  product  as 


UL  =  I®t  (3.57) 

where  I  is  the  N  x  N  identity  matrix  and  i  is  an  L  x  1  index  vector  with  a  1  in  the 
first  position  and  0’s  elsewhere. 

As  an  example,  consider  the  case  where  L  —  3  and  x  e  M4.  The  expansion 
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matrix  U3  G  M12x4  is  given  by 


10  0  0 
0  0  0  0 
0  0  0  0 
0  10  0 
0  0  0  0 
0  0  0  0 

U.3  = 

0  0  10 
0  0  0  0 
0  0  0  0 
0  0  0  1 
0  0  0  0 
0  0  0  0 

which  results  in  y  G  M12. 

Notice  that,  if  L  —  M ,  the  expansion  matrix  is  related  to  the  decimation 
matrix  through  matrix  transposition  as 

Dm  =  (Uif  for  M  —  L\  (3.58) 

in  other  words,  the  operation  of  decimation  and  expansion  are  complementary  pro¬ 
cesses  when  M  —  L. 

3.  Sample  Rate  Conversion  with  Delay 

In  the  analysis  of  various  multirate  implementations,  notably  those  involving 
polyphase  decompositions  (efficient  implementations)  [Ref.  68],  signal  delay  must  be 
incorporated  in  the  system  model  and  forms  another  basic  building  block  of  multirate 
systems  as  shown  in  Figure  3.18.  In  discrete  systems,  a  signal  delay3  is  simply  a  shift 

3Despite  the  name,  there  is  no  restriction  in  considering  advancement  in  sample  index.  This  will 
be  clear  by  the  sign  associated  with  the  delay. 
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of  the  associated  sequence  to  the  right  or  left  by  some  integer  number  of  samples  M 
and  is  represented  by  the  usual 

x[n  —  M]  < — ■»  z~M 

and 

x[n  +  M]  < - ■>  z+M 

where  z~AI  indicates  a  delay  or  shift  to  the  right,  and  z+AI  indicates  a  negative  delay 
or  shift  to  the  left. 

In  conjunction  with  signal  delay,  the  sample  rate  conversion  transformations 
of  (3.50)  and  (3.56)  can  be  generalized.  Given  a  delay  k  G  {0, 1, . . . ,  M  —  1},  then 

y  =  dmx-  (3-59) 

The  constant  matrix  e  M. iVx™  is  called  an  decimation  matrix  with  delay  and  is 
defined  in  terms  of  a  Kronecker  product  as 

D  £>  =  I  ®  $  (3.60) 


where  I  is  the  N  x  N  identity  matrix  and  tk  is  an  M  x  1  index  vector  with  a  1  in  the 
(k  +  l)th  position  and  0’s  elsewhere. 


x[n ] 


♦ 


y[n] 


► 


Figure  3.18.  M-fold  decimator  with  delay. 


In  a  similar  manner, 

y  =  U^x.  (3.61) 

The  constant  matrix  G  ^NLxN  is  called  an  expansion  matrix  with  delay  and  is 
defined  in  terms  of  a  Kronecker  product  as 

=  I  ®  ih  (3.62) 
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where  I  is  the  N  x  N  identity  matrix  and  ik  is  an  Lx  1  index  vector  with  a  1  in  the 
(k  +  l)th  position  and  0’s  elsewhere. 

Consideration  of  (3.60)  and  (3.62)  shows  the  duality  of  the  two  operators,  for 
M  —  L,  and  results  in  the  general  expression 

^m=Pm)T  =  I®*fc  for  /c  G  {0, 1, . . . ,  M  —  1}.  (3.63) 

4.  Linear  Filtering 

It  is  also  useful  to  provide  a  matrix  representation  of  linear  filtering.  Recall 
from  (3.33)  that  a  P  length  causal  FIR  filter  is  fully  specified  by  its  P  filter  coefficients 
(impulse  response).  If  we  denote  these  coefficients  h  G  Rp  as 

hP  =  [  h[ 0]  h[  1]  ...  h[P  -  1]  ]T 

and  the  related  input  data  sequence  x p [ri\  G  Rp  as 

xp[n]  =  [  x[n  —  P  +  1]  ...  x[n  —  1]  x[n]  }T  , 

then  the  output  of  the  filter  at  index  n,  y[n }  G  R  can  be  represented  by 

y[n]  =  hpXp[n]  =  hpXp[n].  (3.64) 


F.  CHAPTER  SUMMARY 

This  chapter  establishes  the  theory  of  multirate  systems  and  provides  the  foun¬ 
dation  upon  which  the  remainder  of  this  work  is  built.  The  concept  of  a  multirate 
system  is  introduced,  with  various  practical  examples.  A  multirate  system  is  formally 
defined,  the  notion  of  rate  is  discussed,  and  the  basic  multirate  operations  of  down- 
sampling  and  upsampling  are  introduced.  A  conceptual  framework  for  the  analysis 
of  multirate  systems  is  developed,  which  enables  a  systematic  extension  of  optimal 
estimation  and  linear  filtering  theory  to  multirate  systems.  The  characterizations  of 
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constituent  multirate  signals  introduced  include  system  rate,  decimation  factor,  sys¬ 
tem  period,  and  maximally-decimated  signal  sets.  Further,  this  chapter  introduces 
the  concept  of  a  system  grid,  which  allows  representation  of  the  various  signals  in  a 
multirate  system  on  a  common  domain. 

The  system  theory  of  multirate  systems  is  also  developed,  including  the  con¬ 
cepts  of  linearity,  shift-invariance,  period  shift-invariance,  and  causality.  In  addition, 
the  input-output  relation  of  a  multirate  system  is  discussed  in  terms  of  the  system  re¬ 
sponse  and  its  associated  Green’s  function  and  is  adapted  to  the  multirate  problem  in 
both  system-level  and  signal-level  representations.  Further,  the  relationship  between 
linear  periodically  time-varying  filters  and  their  linear  time-invariant  equivalents  is 
discussed. 

Finally,  the  basic  multirate  operations  are  analyzed  from  a  linear  algebraic 
point  of  view,  and  matrix  representations  for  the  operations  of  downsampling,  up¬ 
sampling,  sample  rate  conversion  and  linear  filtering  are  presented. 
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IV.  MULTIRATE  OPTIMAL  ESTIMATION 


Signal  or  image  reconstruction  can  be  viewed  as  a  problem  in  signal  estima¬ 
tion,  where  a  related  low-rate  (low-resolution)  signal  or  signals  is  used  to  estimate 
an  underlying  high-rate  (high-resolution)  signal.  From  this  perspective,  the  observa¬ 
tion  signal  or  signals,  and  desired  signal  form  a  multirate  system  and  the  theory  of 
multirate  systems  developed  in  Chapter  III  can  be  used  to  extend  single-rate  signal 
estimation  theory  to  the  multirate  case,  which  is  the  concern  of  this  chapter. 

A.  SIGNAL  ESTIMATION 

“Estimation  is  the  process  of  inferring  the  value  of  a  quantity  of  interest  from 
[typically]  indirect,  inaccurate,  and  uncertain  observations  [Ref.  69].”  A  pictorial 
representation  of  this  concept  is  depicted  in  Figure  4.1,  where  the  source  d  represents 
some  quantity  of  interest  (unknown  parameter,  random  variable,  random  process, 
etc.),  x  are  the  observations,  related  to  d,  and  d(x)  is  the  estimate.  Note  that,  since 
x  and  sometimes  d  are  considered  to  be  random  variables,  d,  which  is  a  function 
of  the  observations,  is  also  a  random  variable.  In  the  field  of  signal  processing, 
these  quantities  of  interest  are  signals,  and  a  major  emphasis  of  research  is  on  signal 
estimation,  where  one  signal  is  estimated  from  some  other  related  signal  or  signals. 
The  desired  signal  may  be  corrupted  by  distortion  or  interference  and  is  usually 
unobservable  (at  least  at  the  moment  when  the  estimate  is  desired).  In  [Ref.  70],  a 
number  of  typical  signal  estimation  applications  are  presented  including:  the  recovery 
of  a  transmitted  signal  from  a  distorted  received  signal,  subject  to  amplitude  and 
phase  distortions  and  additive  white  noise  over  the  communications  channel;  and 
image  restoration  of  an  image  recorded  by  an  imaging  system  that  introduces  blurring, 
nonlinear  geometric  distortions,  and  additive  white  noise. 
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Figure  4.1.  Concept  of  estimation. 


B.  OPTIMAL  FILTERING 

When  the  quantity  of  interest  is  a  random  signal,  optimal  filtering  provides 
a  framework  for  signal  estimation.  Optimal  filtering  is  an  area  of  signal  processing 
study  that  is  concerned  with  the  design  of  filters  to  process  a  class  of  signals  with 
statistically  similar  characteristics  [Ref.  5]  that  are  “best”  in  some  sense,  in  terms  of 
stated  optimality  criteria  (i.e. ,  minimum  mean-squared  error).  The  optimal  filtering 
problem  is  posed  in  the  following  manner.  Suppose  that  a  random  process  x[n]  is 
observed,  which  is  related  to  another  random  process  d[n\  that  cannot  be  observed 
directly.  A  general  expression  for  the  estimate  is 

d[n]  =  0[{x[n]}].  (4.1) 

The  optimal  filtering  problem  is  concerned  with  finding  the  appropriate  functional 
(p[-\  that  provides  the  best  estimate  of  the  desired  signal  d[n\.  When  (p[-\  is  linear, 
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the  functional  is  commonly  referred  to  as  a  (linear)  filter.  This  concept  is  depicted  in 
Figure  4.2. 


Figure  4.2.  General  single-rate  optimal  filtering  problem.  When  </>[■]  is  linear,  the 
functional  is  commonly  referred  to  as  a  linear  filter. 


Typically,  these  filters  are  constrained  to  be  a  linear  function  of  the  observa¬ 
tions  and  hence  can  be  written  in  the  form  (recall  (3.28)) 

OO 

d[n\  —  g[n;  m]x[m],  (4.2) 

m=— oo 

where  the  sequence  g[n\m ],  which  may  be  finite  or  infinite,  is  the  familiar  Green’s 
function.  This  class  of  optimal  filtering  is  called  optimal  linear  filtering. 

Often,  the  optimality  criterion  is  based  on  minimization  of  the  mean-square 
error  (MSE)  between  the  desired  and  estimated  signals,  i.e., 

^2H  =  £(leN2}>  (4-3) 

where 

e[n]  =  d[n\  —  d[n]  (4.4) 

is  the  error  in  estimation  at  the  nth  sample.  When  coupled  with  this  MSE  optimal¬ 
ity  criterion,  optimal  linear  filtering  is  commonly  referred  to  as  Wiener  filtering  in 
recognition  of  the  pioneering  work  of  Norbert  Wiener  [Ref.  71,  72]  in  the  1940’s.1  Of 
great  significance  is  that  only  second-order  statistics  are  required  in  determination 

1Wiener’s  pioneering  work  addressed  the  optimal  filtering  problem  in  continuous  time,  but  is 

easily  translated  to  discrete  time,  as  is  common  now.  Kolmogorov,  in  Russia,  worked  on  the  discrete¬ 

time  problem  for  time  series  [Ref.  73]  and  appears  to  have  preceded  (or  at  least  matched)  Wiener 

in  formulating  and  solving  the  problem  (for  discrete  time). 
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of  these  optimal  filters  [Ref.  61].  Further,  under  the  condition  of  joint  wide-sense 
stationarity  (JWSS)  between  the  input  x[n]  and  the  desired  process  d[n],  the  hlter 
is  shift-invariant  (see  Section  III. D. 2(b)).  For  onr  purposes,  we  develop  the  general 
discrete  FIR  Wiener  filtering  equations,  then  the  filtering  equation  for  JWSS  observa¬ 
tion  and  estimate  signals.  In  order  to  facilitate  such  developments,  the  orthogonality 
•principle  [Ref.  5]  is  stated. 

1.  Orthogonality  Principle 

Theorem  3  (Orthogonality  Principle).  Let  e[n]  =  x[n]  —  x[n)  be  the  estimation 
error.  Then, 

1.  the  optimal  linear  hlter  with  coefficients  h[0],  h[l], . , . ,  h[P  —  1]  minimizes  the 
error  variance  a'l  if  the  hlter  coefficients  are  chosen  such  that  E{e[n\x*[n  —  i\]  = 
£{x[n  —  i]c*[n]}  =  0,  i  —  0, 1, . . . ,  P  —  1,  that  is,  if  the  error  is  orthogonal  to 
the  observations. 

2.  The  minimum  error  variance  is  given  by  <j^m.n  =  £{e[n]d*[n]}  =  £{d[n]£*[n]}. 
The  proof  of  this  result  can  be  found  in  several  places  [Ref.  5,  60,  74,  70]. 

2.  Discrete  Wiener  Filter  Equations 

The  equations  whose  solution  provides  the  optimal  hlter  are  known  as  the 
Wiener-Hopf  equations.  In  developing  the  “ordinary”  (single-channel/single-rate) 
Wiener-Hopf  equations,  let  us  define  the  estimate  in  terms  of  a  linear  FIR  hlter  that 
operates  on  x[n]  and  the  P  —  1  previous  values  of  the  process.  The  estimate  (from 
(3.29))  is  given  by 

p- i 

d[n]  =  h[n\  m]x[n  —  m],  (4.5) 

m— 0 

where  h[n\  m)  is  the  shift- dependent  impulse  response  sequence  of  the  optimal  hlter 
and  P  is  the  associated  hlter  order  (length  of  the  impulse  response  sequence).  If 
x[n  —  i ]  for  i  =  0, 1, . . . ,  P  —  1  represents  any  one  of  the  observations  of  the  sequence 
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x [n] ,  then  by  applying  the  orthogonality  principle  (Theorem  3),  we  can  write 


p- 1 


£{£[n\x*[n  —  *]}  =  £<  |  d[n]  —  h[n\  m]x[n  —  m ]  x*[n  —  i]  >  =  0  (4.6) 


m=0 


or 


p- 1 


^  Rx[n  —  m;  i  —  m]h[n ;  m\  =  P^n;  i 


i  =  0, 1, . . . ,  P  —  1. 


(4.7) 


m= 0 


This  equation  is  the  desired  discrete  Wiener-Hopf  equation  and  can  be  written  in 
matrix  form.  For  example,  for  P  =  3,  we  have 

Rx  [n;  0]  Rx  [n  -  1;  -1]  Rx  [n  -  2;  -2]  h[n;  0]  Rdx [n;  0] 

Rx[n)l]  Rx[n-  1;  0]  Rx[n  -  2;  — 1]  h[n;  1]  =  Prfa:[n;  1]  •  (4.8) 

Rx  [n;  2]  Pz  [n  -  1 ;  1]  Px  [n  -  2 ;  0]  h[n;  2]  [n;  2] 

By  applying  the  second  part  of  the  orthogonality  principle,  the  expression  for 
the  minimum  mean-square  error  can  be  found  as 


p- 1 


ofmiJn]  =  E{e[n}d*[n}}  =  £  <  |  d[n]  —  h[n;  m}x[n  —  m\  d*[n 


(4.9) 


m= 0 


or 

p- 1 

^Ln  M  =  -Rd  [n;  0]  -  /i[n;  m]  [n;  m] .  (4.10) 

m= 0 

If  the  random  processes  are  JWSS,  then  the  statistical  properties  of  the  as¬ 
sociated  system  do  not  change  with  n,  the  estimate  can  be  generated  by  a  linear 
shift-invariant  FIR  filter,  and  the  minimum  mean-square  error  becomes  constant. 

The  matrix  form  of  the  Wiener-Hopf  equation  (4.7)  can  be  derived  directly  if 
(4.5)  is  written  in  vector  form, 


where 


d[n]  =  (x[n])Th, 


x[n 


x[n  —  P  +  1] 
x[n  —  P } 


x[n 


(4.11) 


(4.12) 
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with  x[n]  as  defined  in  Section  II. B. 3,  and 

Mo] 

Mi] 

h\p  - 1] 

or  in  the  general  non-stationary  case 

h[n]  0] 
h[n;  1] 

h[n\  P  —  1] 

Again,  evoking  the  orthogonality  principle,  we  can  write 

R*h  =  rdx,  (4.15) 

where 

R*  =  £{x[n]x*T[n]},  (4.16) 

and 

rdx  =  £{d[n](x[n])*}.  (4.17) 

Equation  (4.15)  has  the  specific  form  (4.8)  (for  P  =  3).  The  minimum  mean-square 

error  can  be  written  as 

<?L in  M  =  Rd  [w;  0]  -  r Ixh* .  (4. 18) 

C.  MULTIRATE  OPTIMAL  FILTERING 

1.  Single-channel,  Multirate  Estimation  Problem 

In  the  single-channel,  multirate  optimal  filtering  problem,  the  observation  sig¬ 
nal  and  the  desired  signal  are  at  different  rates  and  our  goal  is  to  form  an  estimate 


(4.14) 


(4.13) 
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)  x[m\ 

Filter 

d[n] 


Figure  4.3.  General  single-channel,  multirate  optimal  filtering  problem.  Note  that 
the  estimate  and  observation  signals  may  be  at  different  rates. 


of  the  underlying  desired  signal  from  the  observation  signal.  This  notion  is  depicted 
in  Figure  4.3  where  x[m\  represents  the  observation  signal  at  rate  Fx,  d[n]  represents 
the  desired  signal  at  some  other  rate  Fj,  and  d[n\  represents  the  estimate.  A  general 
expression  for  the  estimate  is 

d[n\  =  0[{x[m]}],  (4.19) 

where  the  scalar  estimate  is  a  function  of  the  sequence  x\m\.  Again,  we  seek  the 
appropriate  functional  that  provides  the  best  estimate  of  d[n].  When  </>[•]  is  linear, 
the  functional  is,  again,  referred  to  as  a  linear  filter. 

In  this  work,  we  have  developed  two  approaches  to  formulating  the  required 
estimate.  The  first  approach  is  based  on  a  causal  FIR  Wiener  filtering  model,  the 
second  on  a  non-causal  FIR  Wiener  filtering  model.  The  latter  is  of  particular  interest 
since  it  is  the  basis  for  the  two-dimensional  work  that  follows  (see  Chapter  V). 

a.  Index  Mapping 

The  FIR  linear  filtering  problem  involves  computing  the  linear  combi¬ 
nation  of  some  finite  number  of  observations  to  determine  an  estimate  of  a  desired 
signal  at  a  particular  sample  index,  n  =  uq. 

In  the  single-rate  case,  the  region  of  support  of  the  filter  (see  Sec¬ 
tion  III.D.3)  is  unambiguous.  In  the  case  of  a  causal  filter,  the  filter  uses  x[no\ 
through  a: [no  —  P  +  1]  to  estimate  d[no\-  In  the  multirate  case,  determining  this  re¬ 
gion  of  support  requires  a  way  of  relating  corresponding  sample  indices  between  the 
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desired  and  observation  sequences.  Developing  this  methodology,  referred  to  as  index 
mapping ,  is  the  concern  of  this  section. 

Let  us  motivate  the  discussion  by  considering  Figures  4.4  and  4.5,  where 


d[n ] 
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Figure  4.4.  An  illustration  of  ordinary  causal  FIR  Wiener  filtering  and  the  relationship 
between  samples  of  sequences  d[n\  and  x[n\,  P  —  3. 


the  horizontal  axes  represent  time.  The  first  figure  depicts  single-rate  estimation  of  the 
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Figure  4.5.  An  illustration  of  single-channel,  multirate  causal  FIR  Wiener  filtering 
and  the  relationship  between  samples  of  sequences  d[n]  and  x[m],  P  =  2. 


desired  sequence  d[n]  at  sample  index  n  —  k  with  a  causal  FIR  filter  of  order  P  =  3. 
Since  both  sequences  are  at  the  same  rate,  the  region  of  support,  the  region  encom¬ 
passed  by  the  rectangular  “filter”,  corresponds  to  observations  x[k],  x[k  —  1],  x[k  —  2], 
If  we  consider  the  underlying  discrete-domain  signals,  we  notice  that  both  signals  are 
defined  on  the  same  “signal  domain”,  xI>rd  =  4 >tx-  Since  both  signals  are  defined  on 
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the  same  domain,  the  index  value  k  for  either  sequence  represents  the  same  point  in 
time. 

Figure  4.5  depicts  multirate  estimation  of  the  desired  sequence  d[n]  at 
sample  index  n  —  k  with  a  causal  FIR  filter  of  order  P  =  2.  Here,  we  notice  that  the 
underlying  discrete-domain  signals  are  not  defined  on  the  same  domain.  So,  consider 
a  particular  index  value  k.  In  general,  the  time  at  which  d[k\  occurs  ( kTd )  is  not  the 
same  as  the  time  at  which  x[k\  occurs  ( kTx ).  There  may  be  some  other  index  value 
l  (see  Figure  4.5),  however,  such  that  d[k\  and  x[l]  occur  at  the  same  time  or  very 
nearly  at  the  same  time.  In  an  optimal  filtering  problem,  we  would  generally  like 
d[k\  and  x[l\  to  be  as  close  as  possible.  In  addition,  we  may  also  require  the  filter  to 
be  causal  (see  Definition  23).  Establishing  this  correspondence  is  referred  to  as  the 
index  mapping  problem. 

Let  =  {nTd,  n  G  Z}  represent  the  signal  domain  that  corresponds  to 
an  arbitrary  discrete-domain  signal  with  sampling  interval  Td,  and  dR  =  {mTx\  m  e 
Z}  represent  the  signal  domain  that  corresponds  to  its  associated  observation  signal 
with  sampling  interval  Tx.  Also,  define  a  distance  metric  [Ref.  54]  D[n,  m }  called  the 
index  metric  as 

D[n,  m\  =  | nTd  -  mTx\  =  T\Kdn  -  Kxm\ ,  (4.20) 

where  T  is  the  system  sampling  interval.  The  index  metric  gives  us  a  notion  of 
distance  between  any  two  indices  from  different  signal  domains.  In  fact,  \Kdn—  Kxm\ 
is  the  number  of  system  samples  between  such  indices.  Figure  4.6  illustrates  this 
concept  of  distance  for  two  arbitrary  sample  indices,  no  and  mo-  A  normalized  plot 
of  the  index  metric,  for  Kx  =  3  and  Kd  =  1,  is  depicted  in  Figure  4.7(a).  In  this 
case,  D[n,m]  =  \n  —  3m|,  for  0  <  n,m  <  30.  Each  line  contains  the  locus  of  points 
associated  with  a  particular  value  of  n.  Also,  note  that  the  index  metric  is  identically 
zero  D[n,  m]  =  0  when  n  =  3m. 

The  index  mapping  problem  is  now  stated.  Given  a  discrete-domain 
signal  d  at  rate  Fd  and  its  associated  observation  signal  x  at  a  different  rate  Fx, 
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D(n0,m0)  =  T\Kdn0-Km0 
Figure  4.6.  Notion  of  distance  between  indices  no  and  m q. 


for  a  given  sample  index  n  =  uq ,  associated  with  d,  find  the  observation  sample 
index  m  =  m0,  associated  with  x  such  that  the  related  index  distance  D[no,mo\  is 
minimized,  subject  to  certain  constraints.  Intuitively,  this  makes  sense  since  the  best 
estimate  of  d[n o]  should  typically  result  from  the  closest  observation  samples. 

If  the  problem  of  interest  involves  causal  filtering  (see  Definition  23), 
then  the  required  constraint  on  the  minimization  is 

K 

nT(i  >  mTx  or  n  >  -rr'm.  (4-21) 

Ad 

This  states  that  (on  the  system  level)  the  output  must  not  precede  the  observations. 
Therefore,  the  problem  is  posed  as 


min D(n,m)  subject  ton  > 


mKx 

~1q 


(4.22) 


The  solution  to  this  minimization  is  expressed  as  a  mapping  from  the 
estimation  index  to  the  observation  index  as 


^[■causal  ■  Tl 


nKd 

Kr 


or 


^  3VCcausai(n), 


(4.23) 


(4.24) 
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where  Kd  and  Kx  are  the  associated  decimation  factors. 

This  minimization  is  illustrated  in  Figure  4.7(b),  which  displays  a  plot 
of  the  index  metric,  again  for  Kx  =  3  and  Krj  =  1,  where  the  value  of  n  has  been 
specihed  as  n  =  5.  In  this  case,  the  causality  constraint  requires  that  n  >  3m,  which 
is  indicated  by  the  shaded  region  in  the  lower  figure.  Thus,  the  minimization  pertains 
only  to  that  portion  of  D[5,  m]  contained  within  this  region.  By  inspection,  we  see 
that  when  m  =  1,  the  index  metric  D[5,m]  =  2  is  at  its  minimum. 

Continuing  this  example,  for  Kx  =  3  and  Kfj  =  1,  with  the  causality 
constraint,  we  have  the  following  mappings  from  the  set  of  estimate  signal  indices  to 
the  observation  signal  index  as  shown  in  Table  4.1. 


Table  4.1.  Causal  mapping  from  a  set  of  estimate  signal  indices  to  the  associated 
observation  signal  index. 


If  there  is  no  causality  constraint  (as  in  smoothing  or  image  processing), 
then  the  minimization  is  unconstrained  and  the  problem  is  posed  as 


min  D[n,m  . 


(4.25) 


The  solution  to  this  minimization  is  expressed  as  a  mapping  from  the 
estimate  signal  index  to  the  respective  observation  signal  index  as 


:  n 


(n+m)IU 

Kr 


or 


m  =  Mnc(n), 

where  Kd  and  Kx  are  the  associated  decimation  factors. 


(4.26) 


(4.27) 
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Figure  4.7.  (a)  Normalized  plot  of  D[n,  m]  in  3  dimensions,  (b)  Plot  of  D[n,  m]  versus 
m  for  n  —  5. 
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Again,  from  Figure  4.7(b),  we  can  illustrate  this  minimization,  again  at 
n  =  5,  but  now  with  no  causality  constraint.  In  this  case,  when  the  observation  signal 
index  is  m  =  2,  the  index  metric  is  D[ 5, 1]  =  1,  and  the  index  metric  is  a  minimum. 
Note  that  this  value  of  the  observation  signal  index  m  is  different  from  the  previous 
solution  (a  consequence  of  the  causality  constraint  in  the  former  problem). 

Continuing  this  example,  for  Kx  =  3  and  Kd  =  1,  without  the  causality 
constraint,  we  have  the  following  mappings  from  the  set  of  estimate  signal  indices  to 
the  observation  signal  index  as  shown  in  Table  4.2. 


Table  4.2.  Non-causal  mapping  from  a  set  of  estimate  signal  indices  to  the  associated 
observation  signal  index. 


b.  Single- channel,  Multirate  Wiener- H op f  Equations 

Recall  that  the  input-output  relationship  for  the  general  multirate  fil¬ 
tering  problem  is  given  by  (3.45)  and  is  repeated  here  for  convenience 


y[my\  =  ^  h^[mx]x  ^  ^  ^yl7ly 


lU 


—  m7 


where  /  =  (my) 


My 


This  equation  can  be  written  in  terms  of  the  single-channel,  multirate  estimation 
problem  as 


d[n]  =  ^  [m]x 

m 


~  Kd  n 

-L  Kx  - 


where  l  =  (n)  M  . 


We  can  further  generalize  the  form  of  the  linear  multirate  estimate  if  we  consider  the 
mapping  issues  discussed  in  Section  IV. C.  1(a).  We  can  write  (3.45)  in  terms  of  an 
arbitrary  mapping  function  as 


d[n]  =  V 


m 


(4.28) 
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where  M[-]  is  the  appropriate  mapping  function  (causal,  non-causal) ,  and  l  =  (n) M  . 

With  a  finite  variant  of  (4.28),  in  the  usual  manner,  we  evoke  the 
orthogonality  principle.  Thus,  for  observations  x[M[n]  —  i\  for  i  —  0, 1, . . . ,  P  —  1,  we 
can  write 


fL{e[n]x*  M[n]  —  i  } 

=  £  /  ( d[n ]  —  ^  h^[m]x  M [n]  —  m 


p- i 


m= 0 


for  i  =  0,1,. . .  ,P-1,  (4.29) 


or 

p- 1 

y]  —  m;  i  —  m]h^[m]  =  Rdx[n ;  n  —  M[n]  +i];  %  —  0, 1, . . . ,  P  —  1,  (4.30) 

m= 0 

where  l  =  (n)M  .  This  is  called  the  discrete  single- channel,  multirate  Wiener-Hopf 
equation. 

Applying  the  second  part  of  the  orthogonality  principle,  we  find  the 
following  expression  for  the  minimum  mean-square  error 

p-i 

=  RdW,  o]  -  h[l)  n  -  MN  +  ™\-  (4-31) 

m= 0 

Notice  that  if  the  observation  sequence  is  stationary,  then  the  associated 
correlation  function  is  independent  of  the  mapping  function  and  the  Wiener-Hopf 
equation  can  be  written 

p- 1 

y  Rx[i  —  m]h^[m]  =  Rcix[n\  n  —  M[n]  +  m};  i  —  0, 1, . . . ,  P  —  1,  (4.32) 

m= 0 

where  l  =  {n)M  . 

The  Wiener-Hopf  equation  expressed  in  its  matrix  form,  for  P  =  3, 
when  the  observation  process  is  stationary: 


Rx[ 0]  Rx[- 1]  Rx[- 2] 

h®[  0] 

-  M[n]] 

Rx[  1]  Rx [0]  Rx[—  1] 

h®[  1] 

= 

^x[n;n  -  M[n]  +  1] 

_RX[ 2]  ^[1]  Rx[ 0] 

h«[2] 

Rdx[n;  n  -  M[n]  +  2] 

(4.33) 


where  /  =  (n)M  .  Notice  that  the  correlation  matrix  is  Toeplitz. 

As  in  the  ordinary  Wiener  filtering  problem,  the  Wiener-Hopf  equation 
can  be  derived  in  matrix  form  by  writing  (4.28)  as 


d[n\  =  (x  M[n]  )Th^, 


(4.34) 


where  if 


then 


x  M[n] 


x[M[n]  —  P  +  l] 
x[M[n]  —  P] 


x  M[n] 


x  M[n] 


x[M[n]] 
x[M[n]  —  l] 


5 


x  M[n]  —  P  +  l] 


with  x[M[n]]  as  defined  in  Section  I1.B.3,  and  where  is  defined  as 


(4.35) 


(4.36) 


h« 


h®[0] 

h®[  1] 


h®  [P  -  1] 


(4.37) 


where  l  =  (n)  M  . 

Again,  evoking  the  orthogonality  principle  (Theorem  3),  the  Wiener- 
Hopf  equation  (4.30)  can  be  expressed  in  its  matrix  form  as 


where  l  =  (n)M  ,  and 


Rth(0  =  ? 


dx  3Vt[n]  5 


£{d[n]x*  M[/r]  }. 


(4.38) 


(4.39) 
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Equation  (4.38)  has  the  specific  form  (4.33)  (for  P  =  3).  The  minimum 
mean-square  error  can  be  written  as 


az  \n]  =  Rd\n:  Ol  -?  r 

£min  L  J  “  ’  J 


~T 

dx 


3Vt[r 


(4.40) 


where  l  =  {n)M  . 

c.  Matrix  Approach  to  the  Single-channel,  Multirate 
Wiener-Hopf  Equations 

Another  approach  in  developing  the  multirate  Wiener-Hopf  equations 

involves  the  matrix  representation  concepts  of  Section  III.E  and  results  presented 

in  [Ref.  9].  Again,  consider  a  multirate  system  comprised  of  an  observation  signal 

x  at  rate  Fx  and  a  desired  signal  d  at  rate  Fa,  and  let  the  order  of  the  associated 

filter  be  denoted  P.  Recall  from  (3.10),  (3.14)  and  (3.18)  that  the  system  rate  is 

F  =  lcm (Fx,Fd),  the  associated  decimation  factors  are  Kt  =  F/F, and  the  system 

sample  period  is  IV  =  lcm (Kx,Kd).  The  period  of  the  filters  required  to  form  the 

estimate  d  (see  Section  lll.D.3.b(2),  cyclostationary  period)  is  K,  and  in  the  case  of 

N 

single-channel,  multirate  estimation,  K  =  — — . 

Fd 

For  any  index  n,  the  observation  signal  x ,  represented  by  the  sequence 
x  [n] ,  can  be  expressed 


*(fe) 


k  —  {0, 1, ....  A'  —  1}, 


(4.41) 


where 


x  = 


x[n\ 

x[n  —  1] 
x[n  —  PKX  +  1] 


and  is  a  decimation  matrix  with  delay  (3.60).  In  this  context,  the  decimation 
matrix  with  delay  is  a  mapping  or  transformation  that  extracts  the  samples  in  the 
required  causal  region  of  support  from  the  observation  vector. 


The  optimal  estimate  is  formed  from  a  linear  combination  of  the  ap¬ 
propriate  observation  samples  and  associated  filter  coefficients  and  is  given  by 

dk[n]  =  (xW[n])Th«,  k  =  {0,l,...,K-  1},  (4.42) 

where  the  reversal  operation  x[n]  is  defined  in  Section  II. B. 3.  Note  that  in  this 
formulation,  all  signals  and  computations  are  at  the  system  rate,  and,  as  a  result, 
we  are  solving  a  more  general  problem  since  we  estimate  d,Td  [n]  at  every  point  on  the 
system  grid  not  just  on  the  signal  domain  ^rd-  The  desired  estimate  is  recovered 
by  decimating  the  result  by  K(j. 

In  the  usual  manner,  we  form  the  error  in  estimation  and  evoke  the 
orthogonality  principle.  The  Wiener-Hopf  equation  (4.30)  can  be  expressed  in  its 
matrix  form  as 

Rifc)h(fc)*  =  ?£>*,  k  =  {0, 1, . . . ,  K  -  1},  (4.43) 

or  in  terms  of  single-rate  statistical  parameters  (from  (4.16)  and  (4.17)), 

D^D^hW*  =  D^r*x,  k  =  {0, 1, . . . ,  K  -  1}.  (4.44) 

Applying  the  second  part  of  the  orthogonality  principle  yields  an  ex¬ 
pression  for  the  minimum  mean-square  error, 

a2k  =  i?d[0]  -  f£)Th(fc)*,  k  =  {0, 1, . . . ,  K  -  1},  (4.45) 

where  the  subscript  k  on  the  mean-square  error  reflects  the  periodically  time-varying 
nature  of  the  error. 

2.  Multi-channel,  Multirate  Estimation  Problem 

In  the  multi-channel,  multirate  optimal  filtering  problem,  there  are  multiple 
observation  signals  and  these  signals  are  allowed  to  be  at  rates  different  from  the 
desired  signal.  The  goal  is  to  estimate  this  underlying  signal  from  the  set  of  related 
multirate  observation  signals.  This  notion  is  depicted  in  Figure  4.8  where  there  are 


M  observation  signals  at  rates  denoted  by  m*,  which  may  all  be  different.  A  general 
expression  for  the  estimate  is 

d[n\  =  (f)  [{xi[mi]-,i  =  0, 1, . . . ,  M  -  1}] ,  (4.46) 

where  the  scalar  estimate  is  now  a  function  of  a  set  of  multirate  observation  signals. 
When  the  functional  <p[-\  is  linear,  then  the  estimate  is  formed  from  linear  filtering. 


Figure  4.8.  General  multirate  optimal  filtering  problem  with  M  multirate  observation 
signals. 

a.  Multi-channel  Index  Mapping 

In  this  section,  we  develop  the  mapping  relationship  between  the  esti¬ 
mate  and  observation  sample  indices  in  a  multi-channel,  multirate  system  in  order  to 
determine  the  required  regions  of  support  in  linear  filtering.  In  doing  so,  we  develop 
a  set  of  mapping  functions  for  the  causal  and  noncausal  FIR  filtering  problems.  In 
this  context,  the  index  mapping  problem  is,  for  a  given  index  n,  associated  with  the 
estimate,  to  find  the  M  observation  indices  mll  such  that  the  related  index  distances 
D[n ,  rrij\  are  minimized,  subject  to  certain  constraints,  where  i  —  0, 1, . . . ,  M  —  1. 

If  the  problem  of  interest  involves  causal  filtering,  then  the  required 
constraint  on  the  minimization  is 

K 

nTd  >  m{Ti  or  n  >  — ymj,  for  %  —  0, 1, . . . ,  M  —  1.  (4.47) 

Kd 

Therefore,  the  problem  is  posed  as 

K-i 

min  D(n,  m,i )  subject  to  n  >  —A m *,  for  i  —  0, 1, . . . ,  M  —  1.  (4.48) 

mi  Kd 
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The  solution  to  this  minimization  is  expressed  as  a  mapping  from  the 


estimate  index  to  the  observation  indices  as 


-M-causaZ  •  ^ 


nKd 

K, 


(4.49) 


or 

rrii  =  Mcausai(n),  (4.50) 

where  Kd  and  Kt  are  the  associated  decimation  factors,  and  z  =  0, 1, . . . ,  M  —  1. 

The  multi-channel,  multirate  causal  filtering  problem  is  depicted  in 

Figure  4.9. 

d[n\ 


xj[m]\ 


x2[m2] 


m. 


Figure  4.9.  Concept  of  index  mapping  in  multi-channel,  multirate  FIR  Wiener  filter¬ 
ing. 


If  the  problem  of  interest  involves  noncausal  FIR  filtering,  then  there 
is  no  constraint  on  the  minimization,  and  the  problem  is  posed  as 


min  D[n,  inj  for  i  =  0, 1, . . . ,  M  —  1.  (4-51) 

mi 
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By  examining  the  appropriate  plots  of  D[n,  rrii],  the  solution  to  this 
minimization  is  found  as  a  mapping  from  the  high-rate  index  to  the  respective  low- 
rate  indices  as 

(4- 

or 


Mr, 


:  n 


(n+LfJ)A^ 

K, 


m  =  M  nc(n),  (4.53) 

where  Kd  and  K \  are  the  associated  decimation  factors,  and  i  =  0, 1, . . . ,  M  —  1. 

b.  Multi-channel,  Multirate  FIR  Wiener  Filtering  model 

Recall  that,  for  the  single-channel  FIR  Wiener  filtering  model  (4.28), 
the  estimate  is  written  as 

p- 1 

d[n]  =  ^  h^[mx\x  [M[n]  —  raj  , 

mx= 0 

where  l  =  (n)jv_,  and  M [n]  is  the  required  mapping  function  (causal,  non-causal, 

Kd 

etc.). 

In  the  multi-channel  model,  we  extend  or  generalize  (4.28)  to  include 
multiple  observations  as 

M—l  P-1 

d[n]  =  ^  ^  h!f\m]xi  [Mjn]  —  m] ,  (4.54) 

i= 0  m= 0 

where  p  =  (n)  jv_,  P  is  the  filter  order,  and  Mjn]  is  the  mapping  function  associated 

Kd 

with  the  ith  channel. 


c.  Multi-channel,  Multirate  Wiener- Hop f  Equations 

Using  equation  (4.54),  we  form  the  error  in  estimation  and  evoke  the 
orthogonality  principle.  Thus,  we  can  write 


£{c[n]x*  M i[n]  —  j  } 


M—l  P-1 


£  d[n]  h ^ [s]xr  Mr [n]  —  s\  la:*  M* [ n ]  —  j 


r= 0  s=0 


=  o, 


0  <  i  <  M  -  1,  0  <  j  <  P  -  1  (4.55) 
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or 


M—l  P-1 


y  y  RXrxi\Mr[n]  -  s;  JVC r[n]  -  M i[n]  +  j  -  s]/4p)[s 


r= 0  s=0 

=  RdxAnin  -  Mj [n]  +  j], 

0  <  i  <  M  -  1,  0  <  j  <  P  -  1,  (4.56) 

where  p  =  (n)  jv_ .  These  are  called  the  discrete  multi-channel,  multirate  Wiener-Hopf 
equations. 

As  in  the  ordinary  Wiener  filtering  problem,  the  multi-channel,  multi¬ 
rate  Wiener-Hopf  equations  can  be  derived  in  matrix  form  by  writing  (4.54)  as 

M—l 


d[n\  =  (x,;  [Mi  [n]  ]  )Th^p  , 


i= 0 


where  if 


■[Mi 


n  = 


Xi[Mi[n]  -  P  +  l] 
Xi[Mi[n]  -  P] 


x. 


then 


■[Mi 


n  = 


Xi 

Xi[Mi[n]  -  l] 


Xi[Mi[n\  -  P  +  l] 

with  Xj  [JVC;  W]  as  defined  in  Section  II. B. 3,  and  where  is  defined  as 

hf[  0] 

1] 


h,(P)  = 


h^lP-  1] 


(4.57) 


(4.58) 


(4.59) 


(4.60) 


where  p  =  (n)  n_  . 
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Applying  the  orthogonality  principle  yields: 


M—l 


£{x,;  [Mi[n]]£*[n]}  =  £  <j  (  d*[n]  -  ^  (xr  [Mr [re]]  )*Th{J>')*  j  j  =  0 

0  <  i  <  M  -  1,  (4.61) 


r=0 


and  the  Wiener-Hopf  equations  can  be  expressed  in  their  matrix  form  as 

M—l 


^2  ^  x  TdXi  [Mi  [n]l 

r=0  L  J 


,  0  <  i  <  M  -  1, 


(4.62) 


where  p  =  (n)  n_  ,  and 

Kd 

V  [*[„]]  (4.63) 

Applying  the  second  part  of  the  orthogonality  principle  yields  an  ex¬ 
pression  for  the  minimum  mean-square  error: 


a2  .  [n 

£mm 


M—l 


E,{d[n]e*[n]}  =  £  <|  d[n]  J  d*[n]  —  (xr  [Mr[n]] )*Th^!t 

0  <  i  <  M  -  1,  (4.64) 


r=0 


<y‘l  =  Rd[n;  0]  -  y,  rJXr  [n;  n  -  Mr[n\]  0  <  i  <  M  -  1,  (4.65) 

r=0 

where  p  =  (n)  n  ,  and  its  index  on  the  mean-square  error  reflects  the  periodically 

Kd 

time-varying  nature  of  the  error. 

d.  Matrix  Approach  to  the  Multi-channel,  Multirate 
Wiener-Hopf  Equations 

In  the  multi-channel,  multirate  optimal  filtering  problem,  we  also  de¬ 
velop  a  matrix-based  method  to  develop  the  Wiener-Hopf  equations.  In  this  problem, 
M  observation  signals  are  available  to  form  the  desired  signal  d,  where  M  >  2.  The 
estimate  d  is  formed  by  summing  the  output  of  M  periodic  filters  of  order  P  and  is 
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given  by 


M—l 

dk[n]  =  T:  (x [n\)T h -  ^ ,  0  <  k  <  K  —  1,  (4.66) 

i= 0 

where  K  is  the  period  of  the  filters  required  for  estimation  and  is  given  by  K  =  ( n)  N , 
where  N  is  the  system  period.  For  any  index  n,  the  observation  signal  Xi,  represented 
by  the  sequence  Xi  [n] ,  can  be  expressed  as 

x-fc)[n]  =  D^Xj[n],  0  <  k  <  K,  -  1,  (4.67) 

where 

Xi[n\ 

-  r  ,  Xi[n-1] 

Xj[nJ  = 

Xi[n  -  PKi  +  1] 

with  x, [n]  as  defined  in  Section  II. B. 3,  and  is  a  decimation  matrix  with  delay 
(3.60),  where  Kt  is  the  decimation  factor  associated  with  the  ith  channel,  and  P  is 
the  filter  order.  Again,  note  that,  in  this  formulation,  all  signals  and  computations 
are  at  the  system  rate,  and,  as  a  result,  we  are  solving  a  more  general  problem  since 
we  estimate  d[n]Td  at  every  point  on  4^,  not  just  on  ^Td-  The  desired  estimate  can 
be  recovered  by  decimating  the  result  by  K(j . 

Again,  in  the  usual  manner,  we  form  the  error  in  estimation  and  evoke 
the  orthogonality  principle, 

d[n\  -  j3(xf)[n])Thf)J  |  =  0, 

0  <  i  <  M  —  1,  (4.68) 
and  the  Wiener-Hopf  equation  (4.30)  can  be  expressed  in  its  matrix  form  as 

M—l 

E  R '  =  fS'.  o  <  i  <  M  -  1,  (4,69) 

3=0 
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and  k  =  (n)N,  or  in  terms  of  single-rate  statistical  parameters  (see  (4.16)  and  (4.17)), 


M—l 


E  DigRxixJD^Thf)*  =  Dgr^,  0  <  i  <  M  -  1. 

3=0 


(4.70) 


Applying  the  second  part  of  the  orthogonality  principle  yields  an  ex¬ 
pression  for  the  minimum  mean-square  error  given  by 

M—l 


a „ 


nin [n\  =  £{d[n]e*[n]}  =  £  <j  d[n]  d[n\  -  E  (if')[n])Th-f° 


i= 0 


or 


M—l 


af  =  Rd[n ;  0]  -  ^  0  <  i  <  M  -  1, 

1=0 

or  in  terms  of  single-rate  statistical  parameters, 


M—l 


-1  = 


o]  -  X  hi4’*,  0  <  i  <  M  - 1, 

1=0 


(4.71) 


(4.72) 


(4.73) 


where  the  subscript  k  on  the  mean-square  error  reflects  the  periodically  time-varying 
nature  of  the  error. 


D.  CHAPTER  SUMMARY 

In  this  chapter,  the  problem  of  signal  estimation  is  introduced  from  an  optimal 
filtering  perspective.  Specifically,  we  consider  linear  filtering,  where  the  optimality 
criterion  is  based  on  minimizing  the  mean-square  error  (Wiener  filtering). 

To  begin  the  discussion,  the  Wicner-Hopf  equations  for  single-rate  sequences 
are  reviewed,  establishing  the  foundation  upon  which  the  single-  and  multi-channel, 
multirate  Wiener-Hopf  equations  are  developed.  These  results  are  important  to  the 
methods  of  signal  and  image  reconstruction  discussed  in  Chapter  V. 

Additionally,  the  concept  of  index  mapping  is  proposed,  which  systematically 
describes  the  relationship  between  samples  of  a  signal  at  one  rate  to  those  of  a  signal 
at  another  rate  (indices  in  one  signal  domain  to  those  in  another  signal  domain).  This 
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concept  is  particularly  well-suited  for  developing  the  required  regions  of  support  in 
linear  filtering.  In  detailing  this  relationship,  both  a  causal  and  non-causal  mapping 
transformation  are  developed,  but  the  notion  of  mapping  is  generalized,  providing 
a  basis  for  the  development  of  “generalized”  multi-channel,  multirate  Wiener-Hopf 
equations. 

Finally,  matrix-based  approaches,  using  the  matrix  representations  of  Chap¬ 
ter  III,  are  used  to  develop  the  single-  and  multi-channel,  multirate  Wiener-Hopf 
equations  for  implementation  in  Chapter  V. 
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V.  SUPER-RESOLUTION  SIGNAL  AND 
IMAGE  RECONSTRUCTION 

In  this  chapter,  we  apply  the  multirate  and  optimal  estimation  theory  results 
of  Chapters  III  and  IV  to  the  problem  of  signal  and  image  reconstruction.  We  start 
by  developing  a  1-D  methodology  and  begin  the  analysis  of  this  method  on  a  simple 
test  signal  (a  triangle  wave)  and  evaluate  the  performance  on  this  known  data.  We 
then  apply  the  procedure  to  image  data  when  the  image  is  processed  along  rows. 

We  then  turn  to  the  full  problem  of  image  reconstruction  and  discuss  how 
the  1-D  multirate  theory  and  methods  extend  to  2-D,  present  an  image  reconstruc¬ 
tion  methodology,  and  provide  an  example  of  the  application  of  this  methodology, 
comparing  results  to  other  methods. 

A.  SIGNAL  RECONSTRUCTION 

Consider  the  problem  of  estimating  a  discrete  random  process  d[n] ,  which 
cannot  be  observed  directly,  from  a  set  of  M  related  observation  signals,  represented 
by  {xo [w*o] ,  a;i [mi] , . . . ,  xm-i  [mM-i] } •  These  observation  signals  are  related  to  the 
random  process  d[n]  through  various  forms  of  distortion  and  interference.  Further, 
these  signals  may  be  at  rates  different  from  that  of  d[n] ,  and  observation  signals  at 
the  same  rate  may  also  be  shifted  in  time  with  respect  to  one  another. 

1.  Observation  Model 

In  order  to  facilitate  discussion,  we  present  the  observation  model  depicted  in 
Figure  5.1,  which  represents  the  ith  observation  signal.  In  this  model,  we  consider 
linear  forms  of  distortion.  Notice  that  the  translation  is  along  the  system  grid  (vFjr) 
and  is  indexed  to  the  particular  signal  of  interest  (see  e.g.,  so  and  si  in  Figure  5.2). 
Further,  notice  that  the  downsampling  factor  is  also  indexed  to  a  particular  signal. 
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Distortion 


Translation 


Downsampling 


d[n] 


Additive  Noise 


Figure  5.1.  Observation  model,  where  observation  signals  Xi[mi\  are  derived  from  an 
underlying  signal  d,  subject  to  distortion,  additive  noise,  translation,  and  downsam¬ 
pling. 

If  we  consider  signal  ss[n],  its  associated  downsampling  factor  is  L3,  which  may  be 
different  from  that  of  other  observation  signals. 

2.  Optimal  Estimation 

As  we  have  previously  observed,  if  the  desired  signal  d[n]  and  its  observations 
x,; [mj  are  jointly  wide-sense  stationary,  then  the  linear  filters  required  for  optimal 
mean-square  estimation  are  periodically  time- varying  [Ref.  7].  In  this  case,  we  can 
write  the  estimate  as 

m-  1 

4  [n]  =  Y  (xf }  M  fhf  (5  -  i) 

i= 0 

where  h4  is  a  set  of  time- varying  filter  coefficients  of  length  P?;  and  Xj[n]  a  vect°r 
of  samples  from  the  ith  observation  sequence.  The  periodic  time  variation  is  denoted 
by  the  index  k  where  0  <  k  <  L  —  l,Lis  the  system  periodicity  and  k  =  (n)L  (see 
Figure  5.3).  In  this  analysis,  we  consider  that  the  observations  signals  are  maximally 
decimated  (see  Section  III.C.4)  versions  of  d[n],  where  L*  =  L,  and,  in  this  case,  the 
observation  vectors  can  be  written  as 

x-^[n]  =  D^fc  4)g.[n]  (5.2) 
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Figure  5.2.  Observation  sequences  So  and  si  shifted  by  a  delay  (i  —  0,i  —  1,  respec¬ 
tively). 


and 

s i[n)  =  [si[n],  Si[n  -  1], ... ,  Si[n  -  PiL  +  1]]T  .  (5.3) 

The  decimation  matrix  with  time  delay  extracts  the  appropriate  samples 

from  s i[n]  to  form  each  observation  vector. 

Minimizing  the  mean-square  error  [Ref.  9]  leads  to  a  set  of  Wiener-Hopf 
equations  of  the  form 


p  (fc) 
■^oo 

p  0) 

K0i 

1 

7 

CiT  ^ 

'  o 

1 

tr 

°2 

* 

~(fc)* 

p  (k)*T 
-*R)i 

p  (k) 

''ii 

p  w 

^lL-l 

h[k)* 

= 

~(k)* 

1dl 

p  (fc)*T 
_n'0L-l 

p  (fc)*T 
^lL-l  • 

1 

7 

"  si 

<pc3 

1  ( k )* 

[hL-l\ 

~(k)* 

dL-l_ 

0  <  k  <  L  —  1,  (5.4) 

where  the  time  average  mean-square  error  [Ref.  9]  is  given  by 

°2,=RA  0]-7EEfrhf’-  <5-5) 

k= 0  j= 0 
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Figure  5.3.  Reconstruction  of  the  original  signal  from  an  ensemble  of  subsampled 
signals  based  on  optimal  linear  filtering 


Notice  that  this  is  an  arithmetic  mean  of  the  set  of  0  <  k  <  L  —  l.  The  correlation 
terms  are  defined  as 


=  Df  (5.6) 

=  D^i>£)RiiD^“i>i)*T  (5.7) 

and 

Rd[0]  =  £{d[n]d>]}  (5.8) 

where 

rdi  =  8.{d[n]s*[n]}  (5.9) 

and 

R  ij  =  £{s,;[n]  s*T[n]}.  (5.10) 

Solving  the  multirate  Weiner-Hopf  equations  (5.4)  yields  a  set  of  filter  coefficients, 
which  can  be  used  in  the  estimation  of  d[n]  as  depicted  in  Figure  5.3.  The  application 
of  these  filters  to  the  observation  data  is  illustrated  in  Figure  5.4. 
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Figure  5.4.  Reconstruction  of  the  original  signal  from  an  ensemble  of  subsampled 
signals  based  on  FIR  Weiner  filtering  with  decimation  factor  L  =  3  and  filter  order  P 
=  4.  The  figure  illustrates  the  support  of  the  time-varying  filters  h,-^  at  a  particular 
time,  n  —  15  and  k  =  0  (shaded  circle). 


The  system  equations  can  alternatively  be  derived  using  the  appropriate  map¬ 
ping  function.  In  this  case,  we  find  that  the  causal  mapping  function  can  be  written 


Msr  :  n 


n  —  i 

n 

+ 

k  —  i 

L 

Il\ 

L 

(5.11) 


or 


rrii  =  Msr(n),  (5.12) 

where  L  is  the  decimation  factor,  and  k  =  (n)  L.  Notice  that  the  first  term  I  of 
the  mapping  function  corresponds  to  the  multirate  filtering  scenarios  of  Chapter  IV; 
however,  we  require  an  additional  corrective  term  |_^yj  to  account  for  the  translation 
between  signals. 

Notice  that  if 


n  —  i 

n 

+ 

k  —  i 

L 

Vl\ 

L 

using  the  definition  (2.41)  of  the  common  residue  (n)L,  we  can  write 


n  = 


n 


L  +  (n)L. 
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Then,  by  direct  substitution 


n  —  i 

n  i 

LlJ  L  ( n)L  i 

L 

_L  ~  L_ 

i 

Since  |_^J  is  always  an  integer,  it  can  be  moved  out  of  the  floor  operation  (2.40),  to 
obtain 


n  —  i 

n 

+ 

(n)L 

i 

L 

CL. 

L 

~  L_ 

By  definition,  k  =  {n)L;  therefore,  we  can  write 


n  —  i 

n 

+ 

k  —  i 

L 

CL. 

L 

Example  13.  As  an  example,  for  L  =  3,  we  have  the  following  mappings  from  a 
particular  estimate  index  to  the  corresponding  observation  indices  ( recall  that  for  the 
maximally  decimated  case  with  L  =  3,  i  =  {0, 1,2})  as  shown  in  Table  5.1. 


n 

rrii  =  Msr(n) 

12 

UtJ  +  l¥J}  =  {4>3>3> 

13 

CO 

jjco 

+ 

14 

UtJ  +  l¥J)  =  {4. 4> 4} 

Table  5.1.  Causal  mapping  from  an  estimate  signal  index  to  the  associated  observation 
signal  indices,  for  the  maximally -decimated  case,  L—3. 


These  mapping  transformations  can  be  directly  applied  to  the  multi-channel,  multi¬ 
rate  Wicner-Hopf  equations  (4.62)  and  (4.65). 

3.  Reconstruction  Methodology 

The  process  used  in  the  reconstruction  of  a  signal  from  a  set  of  subsampled 
observations  is  described  as  follows.  First,  a  training  signal,  at  the  desired  rate,  is 
obtained  that  is  representative  of  the  class  of  signals  that  will  be  processed.  From 
this  signal,  a  maximally-decimated  set  of  observations  are  derived  and  the  single- 
rate  statistical  parameters  are  extracted.  Then,  using  the  Wiener-Hopf  equations 
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developed  in  Section  V.A.2,  the  filter  coefficients  are  computed.  With  the  class- 
specific  filter  coefficients,  we  are  able  to  reconstruct  signals  “of  the  same  class”  by 
employing  the  estimate  of  (5.1).  In  other  words,  we  use  training  data  to  develop 
the  required  filter  coefficients  from  a  class-representative  signal  at  the  desired  rate, 
and  apply  these  filters  to  any  set  of  “representative”  observations  at  a  lower  rate  to 
reconstruct  the  desired  signal. 

4.  Application  results 

a.  Reconstruction  of  a  Known  Waveform 

To  evaluate  the  performance  of  the  proposed  method,  two  examples  are 
presented.  In  the  first  example,  a  triangular  waveform  is  considered  for  reconstruction. 
Our  method  was  compared  to  the  method  described  in  [Ref.  31],  which  can  produce 
an  exact  reconstruction  of  the  triangular  waveform  if  the  highest  frequency  terms  are 
left  out.  Both  methods  produce  accurate  results  when  there  is  no  noise  added  to  the 
observation  sequences.  When  a  small  amount  of  noise  is  added  to  the  observation 
sequences,  the  exact  reconstruction  method  fails  to  reliably  reproduce  the  signal  while 
the  method  described  here  continues  to  produce  a  reasonably  good  approximation  to 
the  signal  even  under  severely  noisy  conditions.  The  original  triangular  waveform  is 
shown  in  Figure  5.5  (a).  Also  shown  there  are  the  reconstructed  waveform  (a)  and 
the  mean-square  error  (b). 

The  observation  sequences  are  shown  in  Figure  5.6.  These  were  con¬ 
structed  by  shifting  the  original  sequence,  downsampling  by  a  factor  of  L  —  3,  and 
adding  white  Gaussian  noise.  In  this  particular  example,  a  signal-to-noise  ratio  (SNR) 
of  -4.8dB  was  used.  For  our  purposes,  SNR  is  computed  from  the  ratio  of  signal  power 
to  noise  variance.  Note  that  the  underlying  form  of  the  original  sequence  is  unde¬ 
tectable. 


105 


SNR=  -4.8dB,  Filter  Length  P  =  8,  Decimation  Factor  L  =  3 


Figure  5.5.  Simulation  results  using  optimal  linear  filtering  method  for  reconstruction; 
SNR  =  — 4.8dB,  P  =  8,  and  L  =  3. 


Figure  5.6.  Observation  sequences  of  an  underlying  triangle  waveform  after  being 
subjected  to  additive  white  gaussian  noise  and  subsampled  by  a  factor  of  L  =  3. 
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As  a  precursor  to  consideration  of  the  two-dimensional  image  recon¬ 
struction  problem,  we  applied  the  one-dimensional  methods  described  here  to  image 


data. 

b.  Extension  to  Two-Dimensional  Reconstruction 

We  consider  each  row  of  the  observed  LR  images  as  an  observation  sig¬ 
nal  vector  belonging  to  the  set  {x0[m],  xi[m], . . .  ,XM_i[m]}  (see  Figure  5.7).  Recon¬ 
struction  is  then  accomplished  line-by-line  until  every  row  of  each  image  is  processed. 
In  this  case,  the  original  image  is  depicted  in  the  left  panel  of  Figure  5.8,  and  one  of 


Figure  5.7.  Line-by-line  processing  of  observation  images. 

its  three  subsampled  observation  images  with  additive  white  Gaussian  noise  is  given 
in  the  right  panel.  The  image  depicted  in  the  left  panel  of  Figure  5.9  represents 
the  result  of  applying  standard  nearest-neighbor  interpolation  to  one  of  the  three 
noisy  subsampled  images,  and  the  reconstructed  image  is  shown  on  the  right  panel. 
Although  the  image  is  processed  in  only  one  direction,  there  is  significant  improve¬ 
ment  over  the  interpolated  image.  In  particular,  note  that  edges  of  structures  can  be 
observed  in  many  cases  where  the  interpolated  image  does  not  provide  such  detail. 
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Figure  5.8.  Original  image  (left);  Subsampled  observation  image  L  —  3  with  additive 
white  Gaussian  noise,  OdB  (right). 


Figure  5.9.  Result  of  applying  standard  nearest-neighbor  interpolation  to  one  of  the 
three  noisy  subsampled  images  (left);  Reconstructed  image  (right). 


108 


B.  IMAGE  RECONSTRUCTION 


In  this  section,  we  consider  the  signal  reconstruction  problem  in  two  dimen¬ 
sions,  i.e.,  super-resolution  (SR)  image  reconstruction.  This  analysis  is  an  extension 
of  the  theory  and  application  presented  for  the  1-D  case.  Specifically,  we  consider 
the  problem  of  estimating  a  two-dimensional  discrete  random  process,  which  cannot 
be  observed  directly,  from  a  set  of  related  observations,  at  a  different  sampling  rate. 
Like  the  1-D  problem  posed  earlier,  the  observation  signals  are  related  to  the  random 
process  through  various  forms  of  distortion  and  interference  and  these  signals  may 
also  be  shifted  in  time  with  respect  to  one  another. 

1.  Observation  Model 

In  the  context  of  SR  image  reconstruction,  the  underlying  two-dimensional 
signal  (image)  is  at  a  high-rate  (HR),  and  the  observations  are  at  some  low-rate 
(LR).  The  relationship  between  the  HR  signal  and  a  particular  observation  image 
((i,j) th  —  channel)  is  illustrated  in  the  block  diagram  of  Figure  5.10.  This  observation 
model  shows  that  each  LR  observation  is  acquired  from  the  HR  image  subject  to 
distortion  (typically  blur),  subpixcl  translation,  downsampling,  and  channel  noise. 
We  can  represent  the  observation  model  as 

F«  =  DgFGyDgT  +  Vij.  (5.13) 

The  matrix  F  G  ]^MLixNL2  represents  the  desired  HR  two-dimensional  signal  (image) 
with  (MNL1L2)  pixels.  Its  related  LR  observation  is  represented  by  the  matrix 
F ij  E  MMxAr.  The  matrix  Gij  is  a  linear  operator  that  accounts  for  channel  distortion 
and  matrix  77^  represents  the  channel  noise.  The  parameters  L\  and  L2  represent  the 
horizontal  and  vertical  downsampling  factors,  respectively,  and  the  parameters  i  and 
j  represent  the  horizontal  and  vertical  subpixel  translation,  respectively.  The  matrix 
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Figure  5.10.  Observation  model  relating  the  HR  image  with  an  associated  LR  ob¬ 
servation.  Each  LR  observation  is  acquired  from  the  HR  image  subject  to  distortion 
(typically  blur),  subpixel  translation,  downsampling,  and  channel  noise. 

dl}  is  a  decimation  matrix  with  delay,  defined  in  (3.60),  and  is  used  to  extract  the 
appropriate  pixels  from  the  HR  image  to  form  a  particular  observation  image.  The 
matrix  on  the  left,  incorporates  the  horizontal  downsampling  and  translation  (Li,i) 
while  the  matrix  on  the  right  incorporates  the  vertical  downsampling  and  translation 
(Z/2,  j)-  A  maximally  decimated  set  (see  Section  III.C.4)  of  images  is  obtained  from 

{Ftf},  i  =  0, 1, . . . ,  Li  -  1,  j  =  0, 1, . . . ,  L2  -  1. 


Example  14.  This  example  illustrates  the  operation  of  decimation  matrices  with 
delay  on  a  HR  image.  Consider  Li  =  L2  =  3,  i  =  1  and  j  =  2;  the  HR  image  F  with 
Gi^  =  I  and  rj12  =  0  : 
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The  product  F-ii2  can  be  written 
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The  observation  image  Fij2  is  one  of  nine  (L\xL2)  possible  configurations.  ■ 

2.  Optimal  Estimation 

Consider  the  set  of  LR  observations  { F }  and  the  related  HR  image  F.  We 
desire  to  form  an  estimate  for  the  HR  image  by  some  weighted  sum  of  the  LR  obser¬ 
vations.  The  estimate  can  be  written  as 

F[m,  n2]  =  (fij[ni,n2\,  (5.14) 

i= 0  j= 0 

where  the  expression  on  the  right  represents  the  Frobenius  inner  product  of  the  ma¬ 
trices  ,  defined  in  Definition  12.  Here,  the  matrix  itJ  [n  i ,  n2]  e  is  the  set  of  pixels 

that  form  the  region  of  support  within  the  appropriate  LR  observation,  required  to 
form  this  estimate,  and  matrix  is  the  corresponding  set  of  filter  coefficients.  In 

this  discussion,  this  region  of  support  is  called  the  image  mask  and  its  associated  set 
of  filter  coefficients  is  called  the  filter  mask.  The  Liter  masks  are  chosen  to  minimize 
the  mean-square  error 

£{||F-F||2},  (5.15) 

where  the  norm  is  the  Frobenius  norm;  £{•}  is  the  expectation  operator;  ki  =  (nfi L 
for  0  <  ki  <  Lp  and  (i,j)  represent  subpixel  translation,  with  0  <  i  <  L\  —  1, 

0  <  j  <  L2  —  1,  and  (i,j)  G  Z.  In  the  maximally-decimated  case,  i  and  j  span  the 
entire  set  {0, 1, . . . ,  L\  —  1}  and  {0, 1, . . . ,  L2  —  1},  respectively.  Both  the  set  of  LR 
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image  masks{fjj[ni, n2]}  and  filter  masks  are  further  described  in  Sections 

V.B.2(a)  and  V.B.2(c),  respectively. 

a.  Index  Mapping 

Developing  the  LR  image  masks  {fp  }  involves  mapping  indices  in  the 
HR  sample  index  domain  'hi?  to  those  in  the  LR  sample  index  domain  xVr].  For  a 
given  pixel  intensity  F[rii,  n2],  the  HR  indices  [77-1,712]  map  to  a  set  of  LR  sample 
indices  {[ttii,  which  correspond  to  pixel  intensities  m2]}.  The  indices 

corresponding  to  each  observation  matrix  F,j  are  determined  such  that 

\[ni,n2\  -  [m1,m2]ij\  (5.16) 

is  minimized  for  each  (7,  j). 

This  mapping  can  be  shown  to  be 


[m1,m2\ij  =  [Mi(n1),Mj(n2)], 


(5.17) 


where  M[/r]  is  dehned 


n  —  i  +  [|J 

n 

+ 

k  —  i  +  [|J 

L 

-L- 

L 

mi  =  Msr(n), 


(5.18) 

(5.19) 


where  nTd  e  rriiTi  e  \Eq,  k  =  (n) L,  and  L  is  the  decimation  factor.  Notice  that  the 
first  term  of  the  mapping  function  corresponds  to  the  multirate  filtering  scenarios 
of  Chapter  IV;  however,  we  require  an  additional  corrective  term  to  account 

for  the  translation  between  signals. 

b.  LR  Image  Mask 

The  LR  indices  [mi,m2]jj  for  each  observation  represent  the  centroid 
of  each  of  the  LR  image  masks.  Given  a  desired  mask  size  of  P  x  Q,  each  LR  image 
mask  is  comprised  of  the  P  x  Q  pixels  closest  to  F,:j  [m  1 ,  m2] . 

c.  Filter  Mask 

If  the  desired  HR  image  F  and  its  observations  {F^}  are  jointly  homo¬ 
geneous,  then  the  linear  Liters  H^1,fc2^  required  for  optimal  estimation  are  periodically 
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spatially- varying,  an  extension  of  [Ref.  7] .  This  periodicity  can  be  described  in  terms 
of  the  “phase”  (ki,k2),  where  kt  =  njinodLj.  If  we  define  the  set  of  least  posi¬ 
tive  residues  as  =  0, 1, . . . ,  k  —  1,  then  k\  G  and  k2  G  A l2,  and  all  possible 
combinations  of  phase  can  be  represented  as  A^x  x  A l2. 

Figures  (5.11)  and  (5.12)  depict  the  phase  variation  for  L\  =  L2  =  2. 
In  this  case,  A2  =  {0,1}  and  A2  x  A2  =  {(0,  0),  (0, 1),  (1,  0),  (1, 1)}.  The  spatial 
periodicity  of  the  phase  can  be  observed  by  noting  the  regular  recurrence  of  phase 
terms. 
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Figure  5.11.  Index  representation  to  modulo  representation  with  L\  =  L2  =  2  (note 
the  spatial  phase  periodicity). 


3.  Reconstruction  Methodology 

a.  Least  Squares  Formulation 

In  order  to  determine  the  filter  masks  required  to  estimate 

the  HR  image,  a  least  squares  (LS)  approach  is  employed.  We  identify  the  set  of  all 
HR  pixels  that  correspond  to  a  given  phase  and  denote  this  set  of  HR  pixels  as  the 
matrix  F(fcl,fc2k  This  concept  is  depicted  in  Figure  5.12  where  each  shape  corresponds 
to  a  unique  phase  (circle  (0,0),  square  (0,1),  triangle  (1,0),  and  star  (1,1)).  From 
(5.14),  we  can  see  that 


F[Zi,Z2] 

F[mi,m2] 

F[rai,n2] 


ij 

^(f1J[mi,m2].Hj1,fe2)) 

ij 

5^(fq^i,'«2]-Hq1,A:2)) 


(5.20) 
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Figure  5.12.  Relationship  between  HR  pixels  and  spatially- varying  filter  masks  in 
formulating  the  LS  problem  with  L\  =  L2  =  2. 


where  the  left  hand  expression  is  F(fcl,fe2k  From  this,  we  can  write 

Tp(fcl,fc2)  i®  (JjTjAi’k2) 

—  > 


(5.21) 


where  $  is  the  data  matrix.  This  system  of  equations  is  solved  in  a  least  squares 
sense  for  the  required  set  of  filter  masks  at  each  phase 

b.  Processing  Method 

The  process  used  in  the  SR  image  reconstruction  of  a  set  of  LR  ob¬ 
servations  is  described  as  follows.  First,  a  HR  training  image  is  obtained  that  is 
representative  of  the  class  of  images  that  will  be  processed.  From  this  image,  a 
maximally-decimated  set  of  LR  observations  are  derived  and  then  through  the  least 
squares  methodology  of  Section  V.B.3(a),  filter  coefficients  are  computed.  With  the 
class-specific  filter  coefficients,  we  are  able  to  reconstruct  images  “of  the  same  class” 
by  employing  the  estimate  of  (5.14).  In  other  words,  we  use  training  data  to  develop 


114 


filter  masks  from  a  class-representative  HR  image,  and  then  apply  these  filters  to  any 
set  of  “representative”  LR  observations  to  reconstruct  a  HR  image. 

4.  Application  Results 

In  order  to  evaluate  the  performance  of  this  method  of  image  reconstruction, 
we  process  the  “skyline”  image  depicted  in  Figure  5.13,  subject  to  varying  degrees  of 
AWGN.  The  image  used  for  the  training  process  is  the  204  x  204  pixel  image  segment 
depicted  in  this  figure.  From  this  image,  a  set  of  filter  coefficients  is  derived  that  is 
used  in  SR  image  reconstruction. 


Figure  5.13.  Image  segment  used  to  train  filter. 

The  target  or  object  of  the  reconstruction  is  depicted  in  Figure  5.14.  From 
this  204  x  204  pixel  image  segment,  LR  observations  are  derived  and  are  filtered  using 
the  class-specific  filter  masks.  The  same  level  of  AWGN  is  used  for  the  training  and 
the  target  images. 

Figure  5.15  depicts  three  members  of  the  set  of  LR  observations  with  various 
subpixel  translations.  The  first  image  represents  subpixel  translation  by  one  pixel 
in  the  horizontal  direction  and  no  translation  in  the  vertical  direction.  The  second 
represents  translation  by  one  pixel  in  both  the  vertical  and  horizontal  directions.  The 
third  image  represents  translation  by  two  pixels  in  both  directions. 
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Figure  5.14.  Image  segment  to  be  estimated. 


In  the  remaining  figures,  the  left  panel  depicts  the  SR  image  reconstruction 
using  the  proposed  algorithm  and  the  right  panel  depicts  nearest-neighbor  interpola¬ 
tion  of  one  of  the  LR  observations.  In  every  case,  the  proposed  method  is  superior  to 
the  interpolated  result.  During  these  experiments,  other  interpolation  methods  were 
considered,  including  bilinear  and  bicubic  methods;  again,  the  proposed  method  was 
found  to  be  superior. 

Figure  5.16  compares  reconstructed  and  interpolated  images  for  the  case  of  no 
additive  noise,  downsampling  by  3  in  both  the  vertical  and  horizontal  directions,  and 
with  filter  mask  size  of  3  x  3.  In  this  case,  the  reconstruction  yields  a  result  that  is 
visibly  indiscernible  from  the  target  image. 

Figure  5.17  compares  reconstructed  and  interpolated  images  for  the  case  of  an 
SNR  =  5  dB,  downsampling  by  3  in  both  the  vertical  and  horizontal  directions,  and 
with  filter  mask  size  of  3  x  3.  In  this  case,  we  see  the  effects  of  additive  noise  on  the 
reconstruction.  Despite  some  blurring  of  edges,  details  are  still  discernible. 

Figure  5.18  compares  reconstructed  and  interpolated  images  for  the  case  of  an 
SNR  =  —1.5  dB,  downsampling  by  3  in  both  the  vertical  and  horizontal  directions, 
and  with  filter  mask  size  of  3  x  3.  In  this  case,  the  effects  of  additive  noise  on 
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Figure  5.15.  Downsampled  observation  images  with  subpixcl  translations  (1,  0),  (1, 1), 
and  (2,  2),  respectively;  L\  =  L2  =  3,  P  =  Q  =  3,  and  no  AWGN. 

the  reconstruction  are  quite  deleterious.  Further  blurring  of  edges  is  evident  and 
details  have  become  hard  to  see.  However,  major  features  in  the  image  are  still 
discernible.  Thus,  there  is  a  significant  advantage  in  using  the  proposed  method 
over  interpolation,  where  not  even  major  features  are  discernible.  Note  that  the 
poorer  performance  of  the  interpolation  methods  is  not  unexpected.  In  this  case, 
only  a  single  LR  observation  image  is  used  for  reconstruction  while  in  the  case  of 
the  proposed  method,  multiple,  independent  LR  observations  are  used.  Despite  this 
obvious  disadvantage,  the  interpolation  methods  are  used  throughout  the  literature 
as  a  benchmark  [Ref.  1,  2,  28]  and  we  follow  here. 

C.  CHAPTER  SUMMARY 

In  this  chapter,  the  signal  and  image  reconstruction  problem  is  considered 
from  a  multirate  point  of  view  (Chapter  111),  using  the  multirate  optimal  estimation 
theory  presented  in  Chapter  IV. 

First,  the  problem  of  reconstruction  is  posed  in  one  dimension,  in  terms  of 
a  set  of  observation  sequences  that  are  related  to  a  desired  random  sequence.  An 
observation  model  is  presented  that  models  this  relationship,  accounting  for  linear 
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Figure  5.16.  Comparison  between  a  reconstructed  image  and  interpolated  image; 
Lx  =  L2  =  3,  P  =  Q  =  3,  no  AWGN. 

distortion,  additive  noise,  downsampling,  and  subsample  translation.  Signal  recon¬ 
struction  is  achieved  by  forming  an  estimate  comprised  of  a  linear  combination  of 
samples  of  the  observation  sequences,  and  then  developing  and  solving  a  set  of  re¬ 
lated  linear  equations  for  the  required  filter  coefficients.  These  coefficients  are  used 
to  filter  LR  data.  Finally,  the  results  of  this  proposed  methodology  are  presented  and 
compared  to  other  reconstruction  methods. 

Next,  the  super-resolution  reconstruction  problem  is  considered  and  posed  in 
terms  of  a  set  of  2-D  observation  sequences  that  are  related  to  a  desired  2-D  random 
sequence.  Again,  an  observation  model  is  presented  that  models  this  relationship, 
accounting  for  linear  distortion,  additive  noise,  downsampling,  and  subpixel  transla¬ 
tion.  Image  reconstruction  is  achieved  by  forming  a  linear  estimate  comprised  of  a 
linear  combination  of  pixels  in  each  LR  observation  image  mask  and  developing  and 
solving  a  set  of  related  linear  equations  for  the  required  filter  coefficients.  Formation 
of  these  image  masks  is  based  on  extension  of  the  non-causal  1-D  index  mapping 
results  to  two  dimensions.  The  resultant  filter  masks  are  used  to  filter  LR  image 
data.  Finally,  the  results  of  the  proposed  methodology  are  presented  and  compared 
to  other  reconstruction  methods. 


118 


Figure  5.17.  Comparison  between  a  reconstructed  image  and  interpolated  image; 
L\  =  L2  =  3,  P  =  Q  =  3,  and  SNR  =  5  dB. 


Figure  5.18.  Comparison  between  a  reconstructed  image  and  interpolated  image; 
L\  =  L2  =  3,  P  =  Q  =  3,  and  SNR  =  -1.5  dB. 
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VI.  CONCLUSION  AND  FUTURE  WORK 


A.  SUMMARY 

In  this  dissertation,  new  signal  processing  methods  for  multirate  signals  in  one 
and  two  dimensions  are  developed  and  applied  to  the  problem  of  “super-resolution” 
signal  and  image  reconstruction.  In  super-resolution  processing  a  signal  or  image  with 
high-resolution  is  constructed  by  using  data  from  several  images  of  a  lower  resolution. 

A  significant  contribution  of  this  work  is  the  development  of  the  theory  of  mul¬ 
tirate  systems,  which  provides  the  foundation  upon  which  the  proposed  reconstruction 
methods  are  built.  In  developing  this  theory,  a  conceptual  framework  for  the  anal¬ 
ysis  of  multirate  systems  is  cited,  which  enables  the  extension  of  optimal  estimation 
and  linear  filtering  theory  to  multirate  systems.  Further,  it  leads  to  the  concept  of  a 
system  grid,  which  allows  representation  of  the  various  signals  in  a  multirate  system 
on  a  common  domain. 

System  theory  for  multirate  systems,  also  developed  here  introduces  multirate 
adaptations  of  the  concepts  of  shift-invariance,  periodic  shift-invariance  and  causality. 
In  addition,  the  input-output  relation  of  a  multirate  system  is  discussed  in  terms  of  the 
system  response  and  its  associated  Green’s  function  and  is  adapted  to  the  multirate 
problem  in  both  system-level  and  signal-level  representations.  Finally,  a  method 
is  adopted  to  provide  matrix  representations  for  the  operations  of  downsampling, 
upsampling  and  linear  filtering. 

Another  significant  contibution  of  this  work  is  the  treatment  of  the  problem 
of  signal  reconstruction  from  a  multirate,  optimal  filtering  perspective.  Recognition 
that  signal  or  image  reconstruction  can  be  viewed  as  a  problem  in  signal  estimation, 
where  a  related  low-rate  (low-resolution)  signal  or  set  of  signals  is  used  to  estimate 
an  underlying  high-rate  (high-resolution)  signal,  and  that  the  observation  signal  or 
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signals  and  the  desired  signal  form  a  multirate  system,  motivates  this  work.  In  partic¬ 
ular,  it  suggests  that  extension  of  well-known  results  for  single-rate  systems  could  be 
extended  to  multirate  signals  and  systems.  In  development  of  the  resulting  multirate 
forms  of  the  Wiener-Hopf  equations,  the  concept  of  index  mapping  is  proposed.  Index 
mapping  systematically  describes  the  relationship  between  samples  of  a  signal  at  one 
rate  to  those  of  a  signal  at  another  rate  and  is  particularly  important  for  developing 
the  required  regions  of  support  in  linear  filtering.  In  developing  this  relationship,  both 
causal  and  non-causal  mapping  transformations  are  considered.  The  overall  concept 
is  more  general,  however,  and  provides  a  basis  for  the  development  of  “generalized” 
multi-channel,  multirate  Wiener-Hopf  equations. 

In  applying  the  multirate  and  optimal  estimation  theory  results  to  the  prob¬ 
lem  of  signal  and  image  reconstruction,  a  one-dimensional  method  is  first  explored. 
The  analysis  of  this  method  starts  with  a  simple  test  signal  (a  triangle  wave)  and 
its  performance  is  evaluated  on  this  known  data.  The  procedure  is  then  applied  to 
image  data  when  the  image  is  processed  along  rows.  Finally,  the  full  two-dimensional 
problem  of  image  reconstruction  is  considered  and  an  image  reconstruction  method¬ 
ology  is  developed.  An  example  of  the  application  of  this  method  is  provided  with  a 
comparison  to  results  of  other  methods. 

B.  FUTURE  WORK 

As  fundamental  limits  are  reached  in  future  image  acquisition  systems,  inter¬ 
est  in  signal  processing  solutions  will  continue  to  intensify,  and  in  particular,  the  area 
of  super-resolution  image  reconstruction  will  likely  be  of  great  interest.  During  this 
research,  advances  were  made  in  the  areas  of  multirate  signal  processing  and  optimal 
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estimation  theory,  which  led  to  development  of  the  proposed  reconstruction  methods. 
With  these  hirelings,  other  research  opportunities  became  apparent,  and  are  mentioned 
here. 

In  the  development  of  the  multirate  signal  processing  theory,  several  num¬ 
ber  theoretic  results  were  obtained  to  describe  the  relationship  between  constituent 
signals  in  a  multirate  system.  These  results  were  shown  to  extend  directly  to  the 
second  moment  analysis  of  such  systems,  but  the  second  moment  relationships  were 
never  entirely  understood.  Future  work  in  characterizing  second  moments,  and  even 
higher  moments,  would  provide  further  insight  into  the  behavior  of  multirate  systems. 
Furthermore,  frequency-domain  characterizations  of  such  systems  would  also  provide 
valuable  insight,  and  would  be  a  useful  extension  of  [Ref.  12]  on  bifrequency  and 
bispectrum  maps. 

In  the  analysis  of  performance,  the  2-D  method  developed  here  was  compared 
to  standard  interpolation  methods,  using  a  particular  class  of  images.  In  future 
work,  it  would  be  beneficial  to  compare  the  proposed  method  to  other  “state  of 
the  art”  algorithms  [Ref.  1,  2],  Further,  a  more  robust  analysis  of  the  proposed 
method’s  performance  would  be  achieved  if  several  classes  of  images  were  available 
for  comparison.  Finally,  it  would  be  useful  to  examine  the  effects  of  linear  distortion, 
image  rotation  and  non-Gaussian  noise  on  performance. 

Finally,  in  conducting  this  research,  limitations  were  imposed  on  the  sampling 
rate  of  constituent  signals  of  a  multirate  system.  For  our  purposes,  sampling  rates 
were  constrained  to  be  integer- valued.  The  theory  developed  in  this  work  and  in 
[Ref.  59]  could  be  generalized  to  include  cases  where  real-valued  sampling  rates  are 
considered. 
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